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SYMBOLS 


a = 1,339 = a dimensionless constant 

e l i' e ' > i ~ * w0 vectors 

E(i>) = 47ri»“ • 1/2 4>jj (t^j, V2, ^ 3 ) 

f = longitudinal isotropic correlation coeffic cnt 

g = transverse isotropic correlation coefficient 

HO), HOj, ^ 2 ’ ^ 3 ) = one ‘ ar, d three-dimensional filter functions 

Ky(x) = the modified Bessel function of the second kind of order v 

L = longitudinal integral length scale of turbulence 

Nj, N-> = first and second components of dimensionless spatial frequency 

Ng = dimensionless frequency 

rj = ith component of separation distance 

r = ( fj " + r 2 “ + r 3 “)*/2 

Rjj = cross correlation tensor 

t = time 


T = ambient temperature 
Uj(x,y,z,t) = three component simulated gusts 
Uj(x,y,z,t) = interpolated JAWS winds 
V = aircraft airspeed 

Wj(x,y,z) = three-dimensional block of Monte Carlo simulated turbulence 

x„ „ „ = discrete space function of three variables 
n l' n 2’ n 3 

Xk u ^ = Discrete Fourier Transform (DFT) of x n n 
1 ’ 2’ 3 1> 2» 3 

a(v) = phase of one-dimensional filter function 
P = volumetric expansion coefficient 


vi 



5 jj = Kronecker delta 
6(t) = Dirac delta function 

A^ = lateral gust probe separation in cross-spectra calculations or measurements 

AT = parcel temperature difference from surroundings 

e = turbulent kinetic energy dissipation 

0 = Vt/L = dimensionless time 

©jjfr’l) = one-dimensional spectrum function 

M c ff = effective turbulent viscosity 

r'i = spatial frequencies 

iV = spatial sampling frequencies 
n 

v = 0 |- + v + 1 ~ 

(N^, Ar^/L) = cross-spectra 
^(^1,^2 ^3) = three-dimensional spectrum tensor 
^ijf^l >^2) = two-dimensional spectrum tensor 

p = (r>l“ + ^2")*^" or fluid density. The meaning is obvious from the usage. 


Oj(x,y,z,t) = gust standard deviation 



CHAPTER I. INTRODUCTION 


Approximately every two years in the United States a major wind shear related 
airliner accident occurs killing tens of people. The most recent of these (as of this 
writing) was the July, 1982 crash of Pan Am Flight 759 during takeoff at New Orleans 
International Airport, 156 died aboard the plane and eight otheis were killed on the 
ground as the plane crashed into a subdivision. 

This crash and others like it were caused by wind shear associated with a small 
scale atmospheric phenomenon known as a microburst. lr the past few years two field 
programs were funded to study the microburst. The programs were the Northern 
Illinois Meteorological Research on Downburst (NIMROD) Project and the Joint Airport 
Weather Studies (JAWS) Project. The JAWS Project measured some 70 microburst events 
witn Doppler radar durmg May through August, 1982. Aside from scientific interest, 
several wind shear data sets were subjected to detailed analysis and put into a form for 
use in flight simulator research. These data sets constitute the best wind shear measure- 
ments ever made. 

The JAWS data are presented on a three-dimensional Cartesian grid with grid 
spacings which vary from one case to another but are approximately 200 meters. Hence, 
the JAWS data contain no information on turbulence with length scales shorter than 400 
meters. For nominal landing speeds, some frequencies of importance to aircraft response 
are not contained in the JAWS data. These small scales of turbulence must be added for 
realistic flight simulations. 

During September, 1983 a workshop sponsored jointly by the FAA, NCAR, and 
NASA was held in Boulder, Colorado. The workshop brought together researchers 
directly involved in the JAWS program and potential users of the JAWS data. The users 
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were offered a selected JAWS microburst case and asked what additional information 
they required. From this exchange, several facts emerged. Turbulence was high on their 
list of priorities. Users from private industry look to public agencies such as NASA to 
tell them how to add turbulence to the JAWS data. In the present economic climate, the 
airlines have no money to fund basic turbulence simulation research, nor can they sub- 
sidize the simulator manufacturers to do it. Of the agencies involved, NASA is probably 
best equipped to do the job. 

This document presents an approach to wind simulation which is a significant 
advance in the state-of-the-art. The technique involves the addition of three-dimensional 
Monte-Carlo simulated turbulence to the JAWS data sets. Using this approach, all aero- 
dynamic loads and moments (including roll and yaw) may be calculated from the winds 
simulated over the body of the aircraft. This level of information was previously unavail- 
able from wind simuiaiion models. The spatial model concept, in part, provides the 
answers to the question of how 10 add turbulence to the JAWS data sets. It also serves 
to direct future measurement programs and microcurst research to obt n required turbu- 
lence information. 

In addition, an extension of a previous cross-spectral model based on the von Karman 
turbulence model is presented. Results of the theory are compared with measurements. 
The cross-spectra are a natural part of a three-dimensional simulation. 

A complete explanation of the generation of turbulence, and addition to the 
JAWS data are presented along with FORTRAN computer codes. Background informa- 
tion on turbulence, microbursts, the JAWS data, and Monte-Carlo turbulence simulation 
are also presented. The background material is necessary for a thorough understanding 
of the spatial model. The procedure for generating the turbulence is rather complex and 
an attempt was made to present the material in an intuitive fashion with illustrations and 
geometric interpretations while complex derivations were relegated to the Appendices. 


CHAPTER II. ATMOSPHERIC WIND SHEAR AND TURBULENCE 


Wind shear, or more precisely wind gradients (9u/9x is a wind dilation) have been 
recognized as a cause of aircraft crashes for some time. What was not recognized until 
1976 was the small areal extent of the crash-causing phenomenon. The small scale wind 
shear events were called microbursts by Fujita (11. After the recognition of micro- 
bursts, field programs to study them were done, most notably NIMROD and JAWS. 

One of the achievements of JAWS was the development of microburst data sets for use 
in flight simulation. These data sets are on a relatively coarse grid and therefore do not 
contain information on fine scale turbulence. 

In order to understand the nature of the proposed wind shear model, a descrip- 
tion of microbursts is necessary. The necessary description and the characteristics of the 
JAWS data sets are included in this chapter. These descriptions contain a discussion of 
some microburst models. In anticipation of the need for adding turbulence to the JAWS 
data for the purpose of flight simulation, relevant aspects of turbulence theory are 
presented. 

A. Microbursts 

Fujita defines a downburst as, “a strong downdraft which induces an outburst of 
damaging winds at the surface.” A microburst is defined by Fujita [ 1 ) as a downburst 
of horizontal dimension less than 4 kilometers. A more useful definition is given by 
Wilson and Roberts [2] who define the microburst as a downburst having a differential 
surface velocity greater than 10 meters/sec with the distance over which the velocity 
difference occurs being between 0.4 and 4 km. 

Figure 1 shows the life cycle of a microburst as idealized by Fujita (3]. In the 
figure, descending air meets the ground and spreads out creating strongly diverging 
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surface winds. Microbursts may or may not be associated with rainfall at the surface and 
are accordingly classified as wei or dry. The microburst cases studied during the JAWS 
program show no correlation between surface rain rate and microburst intensity [2]. 

The most intense JAWS microburst had a velocity differential of 48 m/sec (96 kts) and 
was associated with moderate rainfall at the surface [21. The non-correlation of micro- 
burst intensity with rainfall intensity increases their hazard because they can occur in 
apparently benign conditions. 

The final stage of the microburst is called the cushion stage by Fujita. During 
this period surface winds decay and the outflow no longer penetrates to the surface. 

The entire life cycle of the microburst from contact to decay is typically 10 minutes. 
During the mature stage, the maximum wind will occur 50 to 100 meters above the 
surface while the depth of the outflow is typically 600 meters. 

The short duration, high energy, and small length scale of the phenomenon create 
problems for those who would predict, detect, or fly through microbursts. For predic- 
tion, the short duration and random nature of the microburst make forecasting of 
specific events impossible. The best that can b“ done is to forecast conditions conducive 
for the occurrence of microbursts. 

For detection, the small size creates problems. For years at major U.S. airports 
the FAA has operated Low Level Wind Shear Alert Systems (LLWSAS). LLWSAS 
consists of several wind monitoring stations located around the periphery of the airport 
and one center field station. If vector differences of wind velocity between the center 
field detector and any of the other stations exceed a certain level, a wind shear alert is 
issued. Unfortunately, fatal wind shear related crashes have occurred at airports with 
operating LLWSAS, the latest being the crash of Pan Am 759 at New Orleans. One 
problem is that LLWSAS station spacing is larger than the average microburst. The only 
system offering hope of reliable detection in the near future is Doppler radar, but the 
cost of the required system may be prohibitive. Wilson and Roberts [2] have studied 
the requirements of such a system. 
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People encountering microbursts at low levels during takeoff or approach 
frequently find it a once-in-a-lifetime experience. Airliner approach airspeeds are 
typically 75 to 80 m/sec while stall speeds in landing configuration are roughly 60 m/sec. 
With measured airspeed differences of as much as 48 m/sec, microbursts obviously have a 
drastic effect on aircraft performance during takeoff and landing. Figure 2 depicts the 
situation. On approach the aircraft encounters an increasing headwind and begins to rise 
above the glide slope. Pilots have a tendency to throt. back in an attempt to return 
to the glide slope. Passing through the center of the microburst, the plane encounters a 
strong downdraft and an increasing tail wind. A subsequent loss in lift occurs and the 
engines cannot respond quickly enough to save the aircraft. Engines of large air transport 
planes generally require about seven seconds to spool up because of their large inertia. 



Figure 2. Landing aircraft microburst encounter. 


Microburst Models 

The preceding paragraphs contained a brief summary of microburst characteristics 
and effects on aircraft. In the Wind Shear Simulation Workshop in Boulder (Sept. 7-8, 
19S3), simulator users and manufacturers expressed a need for generic microburst models. 
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Because of storage limitations in current simulators, a niicroburst described by algebraic 
equations is needed. We now proceed to a more detailed consideration of inicrobursts 
with discussions of both conceptual models and simple algebraic models. 

Examination of Figure 1 suggests a pair of microburst models. The first is that 
the microburst is a turbulent jet impinging on a surface. The general vertical velocity 
profile through the outflow has much the same appearance as the circular wall jet. 

Because of the transient nature of the microburst, the steady state circular wall jet is 
not a totally satisfactory m^del. In addition, certain turbulence characteristics of a wall 
jet do not seem to agree with measurements [4]. 

A second model requires the occurrence of a source of cooling at cloud base. 

The region of cooling is in effect a momentum source. One mechanism for cooling is the 
evaporation of falling rain. As water droplets evaporate and absorb latent heat, the 
surrounding air is cooled. Figure 1 depicts Fujita’s microburst life cycle. The cool 
parcel interpretation is clear for the contact and outburst stages of the microburst. The 
cushion stage requires some explanation. During the decaying stage of the microburst the 
momentum source at cloud base would be dying. The result is that the falling air would 
be warmer than the air below it so that a stable layer forms near the ground. Because 
the descending air is warmer than air already at the surface, deceleration and spreading 
of the downdraft begin at a higher altitude. 

A third model was devised by Caracena [5] to explain some narrow damage 
swaths observed by Fujita and Wakimoto [6|. Caracena feels that if the jet model is 
accepted a more diffuse damage swath would result. He looked for a mechanism capable 
of maintaining its integrity while transferring momentum over long distances. An obvious 
candidate is the vortex ring. Supporting the vortex ring model are some microburst 
photographs which have a vortex ring appearance. A vortex ring approaching a friction- 
less surface expands and decays similar to observed behavior of microbursts. To illustrate 
his hypothesis Caracena constructed a simple vortex ring generator. If a ring is generated 
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obliquely to a grassy surface, the motion of the grass outlines a narrow “damage swath.” 
Caracena’s arguments concerning jet diffusion are not altogether compelling, however. 

A similar “damage swath” can be created with a hair dryer. 

A question concerning Caracena’s model is the source of the vortex ring which 
requires an impulsive force for generation. One possible source identified by Caracena 
is the collapse of overshooting tops above the anvil of large convective storms. Indeed, 
Fujita and Wakimoto f 6 j have shown some apparent correlation in time between the 
collapse of overshooting tops and the occi.rrence of microbursts. 

Jet, cool parcel (momentum source), and vortex ring models of microbursts are 
conceptual models and full quantitative ramifications of these ideas are not yet developed. 
Two models are developed to the point of algebraic equations, however. These are not 
scientific models, but are intended for use in flight simulators. The equations were 
devised without reference to the equations of flow in one case and with reference to 
continuity only in the other case. 

Bray proposed an extremely simple model for use on Ames Research Center 
simulators. This model is described and critiqued by Elmore [7], The Bray model has 
several shortcomings criticized by Elmore including uniform downdraft source and failure 
to satisfy continuity. Elmore also states that the downdraft diameter is too large for the 
amount of outflow produced. Elmore offers suggestions to make the Bray model more 
realistic. These suggestions include modification of the source terms and alteration to 
satisfy continuity 

A more sophisticated model was proposed by Zhu and Etkin, 1983 [8] which is 
an ideal fluid model. It involves the flow induced by a doublet disc of variable strength 
together with its image. With this model, Zhu and Etkin were able to demonstrate the 
excitation of aircraft phugoid frequencies by microbursts. 

The Bray model, though admittedly crude, does provide a reasonable simulation. 
With the addition of Elmore’s suggestions, simulations should be improved. Of the three 
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models described, Zhu’s and Etkin’s model is probably the best. Neither the Bray nor 
the Elmore model were intended as scientifically accurate pictures of microbursts and as 
a result some of the physics not considered essential for simulation are neglected. 

Neither of the three show a microburst front and consequent abrupt shear at the leading 
edge. The equations which can be made functions of time as well as space are assumed 
separable, i.e., f(x.t) = X(x) T(t). This is not a realistic picture since it implies that the 
microburst intensifies and decays uniformly at every point. 

The preceding discussions have dealt with conceptual and algebraic pictures of 
microbursts. In the next subsection, characteristics of real JAWS measurements are 
described. 

JAWS Data Sets 


At the present time three of the JAWS data cases have undergone multiple 
Doppler analysis to determine three component winds. These three cases are summarized 
in Table 1. 

Table 1. JAWS Dual Doppler Cases 


Date 

Ax' (m) 

Ay“ (m) 

June 29, 1982 

300 

300 

July 14, 1982 

200 

200 

August 5, 1982 

150 

150 


1 . Ax is the east-west grid spacing. 

2. Ay is the north-south grid spacing. 

3. Az is the vertical grid spacing. 


Az 1 2 3 (m) 

AV max < m ' /sec) 

Number of 
x x y x z 
grid points 

250 

25 

60 x 60 x 9 

150 

30 

60 x 60 x 1 1 

250 

30 

81x81x9 
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The August 5 case was selected for release to selected companies because of its 
high resolution and intensity. The data were summarized in a document [9] which con- 
tains a plot of low level winds showing a sinuous jet emitted from the northeast quadrant 
of the microburst. This feature also appears in the winds in the next level up but 
gradually disappears at higher levels. The jet seems to be a real phenomenon and if so 
it raises some questions. The main question is, “Why is it there?” Is this a reflection of 
a convoluted source up at cloud base? This idea seems highly unlikely since diffusion 
tends to reduce source convolution in nature. The alternative is that something is 
happening at the surface. In other words the jet may be a result of topography or of 
some surface roughness anomaly. Reference to a Denver map shows that the jet occurs 
in the near vicinity of the South Platte River basin and the jet flow is in the same 
general direction as the river. The implication is that the cold stable air is being 
channeled along the river basin. 

The JAWS data is generally displaced freely relative to airports in simulations. 

If the jet is a surface effect, freely moving the data to flat areas surrounding airports 
may not be a reasonable procedure. 

Figure 3 is a vertical section through the July 14 case which is a typical asym- 
metric microburst embedded in a wind field. The center of the section shows a strong 
downdraft. On the right side of the downdraft is a slight updraft. The width of the 
downdraft is roughly 1 km. This downdraft is just off center of a moderate to heavy 
rainshaft and this microburst would be classified as wet. 

Microbursts weren't identified until 197b and many questions concerning their 
nature remain unanswered. What is the best model of a microburst (jet, momentum 
source, vortex ring, etc.). What is the effect of source size and strength on outflow 
intensity? What effect does source height have on the microburst? What effect does 
surface topography and/or roughness have on the outflow pattern and damage swath? 
What degree of kinematic organization exists in a microburst (vortex ring)? How do 
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turbulence characteristics (length scale and intensity) vary through a microburst? This 
latter question is crucial for turbulence simulation. JAWS has provided partial answers 
to some of the above questions. Most require additional measurements preferably in the 
controlled environment of a laboratory. Appendix A is devoted to laboratory approaches 
for obtaining hese answers. 


B. Isotropic Turbulence 


Atmospheric turbulence is neither isotropic, homogeneous, nor Gaussian. Despite 
this fact at the higher frequencies and shorter length scales, turbulence tends toward 
isotropy. Isotropic turbulence is the simplest form of turbulence. In the discussion of 
the spatial wind model, a means for separating the longer wavelengths and non- 
homogeneities will be discussed. The need for the separation arises from the difficulty 
of generating nonisotropic turbulence with Monte-Carlo techniques. In anticipation of 
a later need, some useful properties of isotropic turbulence are presented in this subsec- 
tion. 


The description of the kinematics of homogeneous isotropic turbulence begins 
with the definition of two functions, the longitudinal and transverse correlation 
functions. These two functions, defined in Figure 4. are functions of separation r, and 
time. A similarity form which is a function of r/L only is assumed. Time dependence is 
contained in L = L(t). The best known mode's of f and g are the von Karman and 
Dryden models, both of which can be written as functions of r/L only. By virtue of 
continuity a relation exists between f and g. 
r df 

8 = f+? - a) 


Different investigators have used different functions to fit the measured correla- 
tions. The two most famous correlation models are by Dryden and von Karman. The 
Dryden correlations are given by 
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Figure 4. Longitudinal and transverse correlation functions. 


f(r) = e~ r / L 

(2) 

g(r) = (1 - r/2L) e~ r / L 

(3) 


where L is the longitudinal length scale of turbulence ( / f(r) dr). These equations 

J Q 

have the advantages of simplicity and rational spectra. Simplicity is an obvious advan- 
tage, but the advantages of rational spectra are related to mechanisms for simulating 
turbulence by Monte-Carlo methods. The main disadvantage of the Dryden model is that 
it does not Fit measured data quite as well as the von Karman model. 

The correlation functions of the von Karman model are given by the following 
equations'. 



where 

T(x) = gamma function 
a = a constant = 1.339 


I 


Kj/x) = the modified Bessel function of ihe second kind of order v . 
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The correlation functions are important because the complete three-dimensional 
correlation tensor can be expressed in terms of these two functions. 



(f-g) 




where 


5 -j = 1 i = j Kronecker Delta 
= 0 i ¥=} 


( 6 ) 


The above result was obtained by direct calculation by von Karman and Howarth [10]. 
The tensor function can be used to calculate correlations between velocities oriented 
along any two unit vectors, say ej j and e-»j. The scalar correlation between these two 
velocities is given by 


R 12 - R ij e l,i e 2j (7 

In the above equation the Einstein summation convention is in effect. 

Figure 5 gives a geometric interpretation of the preceding discussion. Since Ry 
is an isotropic tensor, R^ is invariant with respect to rotation, i.e., if the three vectors 
v j , vi, and r are rotated in any manner through the isotropic, homogeneous turbulence 
so that they maintain the same orientation with respect to each other (rigid body rota- 
tion) then Rj2 ‘ s unchanged. 

The correlation tensor has an equivalent representation in frequency space. 


oo oo oo 


♦jjC*'! Wi* = Iff R ij( r l’ r 2< r 3) e j27r(r l^ +r 2' ; 2 +r 3^3) drj dr2 df3 


-OO — oo —oo 


!t 


E(») 


A ^ 


- (v^ j/.) 

1J 


( 8 ) 
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where pj, i> 2 , and ^3 are the three spatial frequency components 

-) -> 7 1/2 

y = j “ + ^2" + ^3 ) 

and EO) = 4m>" • 1/2 d’yfr'j, ^ 7 , ^ 3 ) and the Einstein simulation convention applies. 

r 3 



Figure 5. Geometry of the cross-correlation tensor. 

As in the case of correlation, <hjj is a tensor function and a scalar function is 
obtained as before. 


W 1^2^3> = e l,i e 2 J (S 

In the above expression ej j and C 2 j are the same unit vectors which correspond 
to vj and V 2 in Figure 5, and 4 >j 2 is the scalar cross-spectrum corresponding to vj 
and vj- 
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EO) has representations corresponding to the Dryden and von Karman models. 
For the Dryden model the equation is 


E(t>) = 16 a“ L 


(27rLi’) 4 


[1 + (2ttL^) 2 ] 3 


( 10 ) 


The corresponding von Karman equation is 


1 10 T 

E(y) = a ~ L 


(27raL^) 4 


[1 + (27raL^)^" 1 


^ ' 17/6 


( 11 ) 


From these two equations <f>jj for the two models can be calculated. For the Dryden 
model the result for the auto-spectra is 


4> i i (*'l’ l; 2' l '3) = 


64 o~ L 3 7 r 3 {v 3 - vf) 
[1 + (2?rLn) 2 ) 3 


( 12 ) 


In this equation and the next, the Einstein summation convention is suspended. 
The von Karman model is 


*ti v \’ v 2 ^ 3 ) = 


440 n 3 ° 2 a4 l5 {v2 ' "i 2 ) 
9 [1 + (27raD') 2 J 17/6 


(13) 


\ 


The spectra for the vertical velocity component, equations (12) and (13), can be 
rewritten for the Dryden and von Karman models, respectively, as follows: 


Dryden Model 

^33(^1 ,^ 2 ^ 3 ) 


64 0 “ L" 1 7i (r>j 2 + ^ 2 ) 

[1 + (27rL) 2 (r>j 2 + r> 2 2 + J^ 2 )] 3 


(14) 
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von Karman model 


. , _ 440 7r^ "’4.5 

®33 {v \’ v 2' v 3^ 5 a L 


(t’j ~ + t’2 ) 


[ 1 + (27raL)“ (v\ - + i i>o“ + ^ 3 “)! 1 7 ^ 6 


• ( 15 ) 


Equations (14) and (15) show rotational symmetry about the axis, i.e., they 
are functions of p“ = nj~ + u^~ and Working with the von Karman spectrum. 
Equation (15), we can solve for = ^3(^33 , p). The result is 



where Cj and C-> are functions of a and L. 

Equation (16) is an equation for a surface of constant in wavenumber 
space. Values of Cj and C-, for the corresponding spatial frequency spectrum are 
developed in Appendix B. 

Equation (16) is the equation of a toroid. C oss sections of it ate plotted in 
Figure 6 and a perspective of the surface is depicted in Figure 7. Surfaces of constant 
$33 are found if the curves of Figure 6 are rotated about the v 3 axis. A result similar 
to (16) can be obtained for j and with rotational symmetry about the v j and v~> 
axes, respectively. 

Figure 6 seems to imply that $33 is singular, i.e., it approaches infinity near 
the origin. In fact <1*33 does have a finite maximum value. 

Three-dimensional correlations and cross-spectral tensors are extremely difficult 
to measure. Actually, one-dimensional correlations are measured and corresponding 
one-dimensional spectra calculated from the correlations, i.e. 


R 


j(r j ,0,0) = / 


OO OO 

If 4 > ij(t'], I '2’ I '3) d "3 


— j27rt»i r j 
e 1 1 dt^i 





tMMIMII Mu a>* 


(i)‘ 


, x 64 t r 9 <1 

^33(^1 ," 2 ) = — a(aLr 
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(^1“ + v 2> 

9 " 9 9 tTT ( 21 ) 

[1 + (2ffaL)" (t’j' + ^2^)1 7 ^ 


©j ](*",)= 2a“L 


[ 1 + (27raLt'j)^]^ 


O ^ 

1 + - (27raL^j)^ 

0 2 2(^,) = ©33(^1) = a 2 L — — (23 

(1 + ( 2 ™^, ) 2 ] 11/6 

The results from this subsection are used in the next subsection to - alculate a 
theoretical cross-spectral model. In Chapter VI results from the preceding d^-ussion are 
used to calculate filler functions for generating three-dimensional turbulence. 


C. A Theoretical Cross-Spectral Model and Comparison with Measurements 


Elements of the previous discussion concerning isotropic turbulence will now be 
used to develop a cross-spectral model of turbulence. The assumption of frozen, von 
Karman turbulence is made. Previously Houbolt and Sen [12] did a similar analysis to 
obtain the cross-spectrum corresponding to the vertical velocity component. An equiva- 
lent approach is used to extend their analysis to the longitudinal and lateral cross-spectra. 

In the analysis that follows it is convenient to non-dimensionalize the spectra. 
Without non-dimensionalization. annoying conversion constants must be carried along so 
that the units are correct. Equations (19) through (21) are non-dimensionalized using 
the transformations Nj = v j L. 


^1^9) = ^ 




Applying this transformation to (19) through (21) yields: 

1 + ^rraNj) 2 + ^ (2traN 2 ) 2 

*1 i< N i’ N 2> = j (oa > 2 9 — 9 TTh 

3 [1 + (2war (N, 2 + N 2 2 )] 7/3 
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The frequency transformation, Nj = ►'jL implies a transformation in space, X = 
rj/L, Y = r->/L. Then 


oo oo 

f r +}2ir(N|X+N«iY) 

Ry(X>Y) = j J* ij(N 1 .M 2 ) e 1 - dNjdN 2 


(28) 
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The object we are trying to achieve is the cross-correlation between the velocity at one 
aircraft wing tip at one moment and the velocity at the other wingtip at a later moment. 
Assuming frozen turbulence. Figure 8 applies. Making the transformation 0 = Vt/L gives 

N-,) e +j 27r ( N l 0+N 2 Ar 2/ L > d Nj d N 2 . (29) 

Equation (29) is an expression for the desired cross-correlation. To obtain the desired 
cross-spectra, the Fourier transform of (29) is taken. 
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The exponential integral with respect to 0 above is the Dirac delta function fifNj-N^), 
so the final result is 
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Figure 8 . Illustration of cross-correlation distance. 
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Using Equation (30) along with Equations (25) through (27) gives the final expressions 
for the three cross-spectra. 
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In these expressions the evenness of the integrand was used. The above integrals can be 
evaluated using Filon’s method [11]. The results are plotted in figures 9 through 1 1. 

Comparison is made with data collected during the NASA B-57B Gust Gradient 
Program. The B-57B is instrumented with three component gust probes on each wing 
tip and at the nose. Figures 12 through 14 show comparison with data collected during 
a horizontal flight. Figures 1 5 through 1 7 depict a comparison for flight on a simulated 
ILS approach (three degree glide slope). Both measurements were taken within the 
planetary boundary layer where turbulence is not really isotropic. 

In some of the figures, the measured spectra fall between the expected Ar->/L 
curve and the Ar^/L = 0 curve. This result is quite apparent in figure 14, for example. 

A possible explanation concerns the assumption of the stationarity of the correlation 
function. In the real atmosphere, the correlation falls off with time as well as with 
space. This increased fall off in space and time corresponds to a slower falloff in the 
frequency domain as observed in the figure. 

Some of the measured data show an unexpected foot (flattening out) at the 
higher frequencies. Possible sources of this error include aliasing and the addition of the 
white noise to corresponding velocity components measured at each wing tip. The 
correlations would show a Dirac delta function spike which in the frequency domain 
corresponds to a small but finite constant. When plotted on log paper the spectra would 
show a characteristic foot. Simultaneous white noise could be introduced by terms from 
the INS system used to remove the aircraft motions from the data. 

The importance of a cross-spectral model is that moment spectra can be 
calculated from the cross-spectra [13], 
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Figure 14. Comparison of measured and computed w-component 
cross-spectra for a level flight case. 
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Figure 15. Comparison of measured and computed u-component 
cross-spectra for the simulated ILS approach case. 
















CHAPTER III. REVIEW OF MONTE CARLO TURBULENCE SIMULATION 


For the purposes of this study, turbulence simulation will mean Monte Carlo 
turbulenee simulation Monte Carlo turbulence simulation should not be confused with 
efforts to close me time averaged equations of motion. Monte Carlo turbulence simula- 
tion is defined as a procedure whereby a noise process is filtered by a linear or nonlinear, 
analog or digital filter to obtain an output with certain of the statistical properties of 
turbulence. While Monte Carlo simulation and numerical simulation of turbulence have 
some things in common, the intent of each is quite different In the case of the former, 
mean velocity profiles are known and turbulence is to be added to create a realistic How 
field. In the latter case mean velocity profiles are to be calculated. Monte Carlo simula- 
tion strives to put small scale perturbations back into the flow field while numerical 
simulation attempts to remove the small scale details. In this sense the two are opposite 
operations. 

Monte Carlo turbulence simulation apparently began in the mid-fifties. The basic 
idea, which is employed in the present study, is depicted in Figure 18. Gaussian white 
noise, with Dirac delta function autocorrelation and corresponding constant spectral 
density is input to a linear filter. The linear filter has a transfer function H(u) which is 
selected to give the desired output spectrum. Because the filter is linear, and the input 
Gaussian, the output is Gaussian. Using the notation of Figure 18, the output spectrum is 

QJp) = IHOOI 2 © n 00 (34) 

In (34), © w (^) is the known, desired form of the output spectrum, and © w (*0 is a con- 
stant for all frequencies. H^) is the unknown with modulus 
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Figure 18. Simple Monte Carlo turbulence 
simulation. 


From Equation (34) the realization ot H(r) which gives the desired 
output spectrum is nonunique. In fact, H(i’) has the general form 

H(v) = |HU0I e + j“ (z '> . (35) 

where u(i’) is the phase of H(t'). H(r>) has an infinite number of realizations correspond- 
ing to an infinite number of choices for the phase. One choice of phase is a(r') = 0 for 
all i;. The choice of the phase docs not effect the output spectrum. 

Choice of H(r>) is complicated by the fact that no noise source, neither analog 
nor digital, is completely \ hite. In the digital case, which is the concern of this report, 
the designated sampling frequency determines the cutoff frequency of the noise. The 
nonwhiteness of the noise source results in a factor wliich is multiplied by the expression 
for Hfr 1 ). This aspect of the simulation is discussed in more detail in Chapter VI. 

Not all turbulence simulations involve the us^ of linear filters. Reeves and nis 
colleagues in a series of papers [14-16] developed a nonlinear filter model. Figure 19 
depicts the block diagram of one version of the nonlinear filter used by Reeves and 
associates. In the figure, the parameter R can be changed to modify the kurtosis of the 
output turbulence. In effect R changes the patchiness of the simulated turbulence. The 
functional form* of the three filters H a , H^, and H . are chosen to give the desired output 
spectrum. 
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Figure 19. Reeves' non-Gaussian simulation model. 


In general, the output form of the probability density eannot be determined for 
nonlinear Filters. Reeves judiciously constructed his Filter so that the output probability 
density function could be calculated. 

Another advance in turbulence simulation arose from the need to include inter- 
level coherence m simulation. Fichtl and Perlmutter [17] developed a multi-filter system 
to incorporate interlevel coherence. In this model a series of white noise sources is 
filtered, added together, and the result filtered to obtain the ti rbulence. A schematic of 
the method is given in Figure 20. Each of the Filters Dp through Dp are effectively 
height dependent phase shifts. The desired coherence is obtained exactly as the number 
of Filters and noise sources approaches infinity. In practice, the desired coherence can be 
approximated to any desired accuracy by choosing p large enough. The filter H in Figure 
20 gives the desired output spectrum. Since the block diagram of Figure 20 represents a 
linear Filter and since the inputs are Gaussian, the output will also be Gaussian. 

For the purposes of this study, one of the more significant advances in turbulence 
simulation technology was made by Fichtl [18J. In this report, Fichtl generated non- 
dimensional turbulence. By this approach, one nondimensional turbulence record could 
be generated, its spectrum checked, and then be used for simulation of any flight profile. 
Fichtl’s method was based on the Dryden spectral model. A similar approach was used 
by Tatom, et. al. [19] for the von Karman model. 














36 



Figure 20. Interlevel coherence model [17], 


The preceding discussion was a brief history of some turbulence simulation 
approaches which have some relevance to the current model. For a broader discussion 
of Monte Carlo turbulence simulation, the interested reader is referred to three papers 
by Dutton, et. al. [20], Fichtl, et. al. [21], and Wang and Frost [22], The report by 
Dutton, et. al. gives a good historical perspective of turbulence simulation and presents 
an unusual nonlinear approach. The paper by Fichtl, et. al is a general survey of simu- 
lation methods up to 1977. The paper by Wang and Frost presents codes and discussions 
for several of the methods discussed previously. 

Before leaving the area of Monte Carlo turbulence simulation, a brief discussion 
relevant to the present model is presented. In Chapter II, one-dimensional spectra for 
the von Karman and Dryden models were given. Real turbulence shows a rolloff 
as predicted by Kolmogorov and the von Karman model, however, the von Karman 
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model lacks the convenience of a rational O"^ re' t .) Dryden model. For rational 
spectral models, the simulated turbulence can be generated with a difference equation 
of the form 

w i+l = W w i> w i-l> • • ■ • n r n i-l> • ) (36) 

Equation (37) is a digital simulation model where the i+lst turbulence point generated 
is a function of the previous turbulence points and noise source points. For the Dryden 
model (37) can be written as 

w i+j = c i w i + Co w i_i + dj n i + dT n i_i , (37) 

where Cj, C->, dj, and d^ are parameters depending on the sampling rate, airspeed, 
turbulent intensity and length scale of turbulence. 

Wang and Frost [22] describe a simple rational model which approximates the 
more realistic von Karman spectra. This approach can be carried to any level of com- 
plexity as long as care is taken so that all poles lie within the unit circle in the Z- 
transform plane. This stipulation assures system stability. 

When a recursion relation such as Equations (36) or (37) can be derived, it is much 
superior computationally to other available methods. Otherwise the turbulence must be cal- 
culated either with a convolution or with a Discrete Fourier Transform (DFT). When turbu- 
lence in two or more dimensions is generated, the DFT or convolution must be used because 
Z transform theory is not well developed in more than one dimension. The only time 
it can be used is when ihe function to be transformed is separable so that the two- 
dimensional Z transform reduces to the product of two one-dimensional transforms. 
Unfortunately, spectral models are functions of v = {vy + + v y) 1 which is not 

separable. 

In this chapter a brief review of Monte Carlo turbulence simulations was pre- 
sented. The review was by no means exhaustive, but was confined to those approaches 
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which are relevant to the current method. With the exception of the paper by Fichtl 
and Perlmutter [17] none of the reported studies in any way accounted for variations 
of turbulence in more than one dimension. Fichtl, et. al. |21] discuss in passing the 
generation of two-dimensional turbulence. For truly realistic turbulence simulation, 
variation of the turbulence in all three dimensions is necessary. The generation of 
three-dimensional turbulence is the subject of Chapter VI. 



CHAPTER IV. DESIRABLE FEATURES OF A SIMULATION OF 
ATMOSPHERIC WINDS AND TURBULENCE 


Monte Carlo simulation models of atmospheric winds and turbulence have several 
desirable features. Some of these are listed and discussed in the following paragraphs. 

1. Realism: The model should be realistic in the sense that turbulence generated 
by Monte Carlo methods should have as many of the statistical and spectral characteris- 
tics of real turbulence as possible. The turbulence should have the right “feel.” A 
common complaint of pilots who do some of their training on simulators is that the 
generated turbulence is not realistic. This absence of the correct feel has been attributed 
to failure to simulate turbulent intermittancy or patchiness [14], Another factor may be 
the failure to simulate variation of winds across the body of the aircraft. By this failure, 
the main contributor to roll and yaw moments are neglected. Six degree of freedom 
motion simulators (controlled by hydraulic actuators) exist so that roll, pitch, and yaw 
motions could be simulated if corresponding moment information were available. 

Finally, a simulation model should contain all frequencies which cause a significant 
response in the man-aircraft system. 

2. Computational Speed: The wind and turbulence model should be computa- 
tionally efficient so that real time flight simulations are possible. A computer working 
with a flight simulator must process incoming analog signals, convert digital output 
signals to analog response, calculate turbulence and aerodynamic loads all in real time. 
Many operational flight simulators are operated by minicomputers and computational 
speed is a must for pilot-in-the-loop simulations. Research simulators may have one or 
more high speed state-of-the-art computers to run them but even these may not be able 
to handle a real time turbulence simulation if the model is too complex. 
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3. Flexibility: Several types of flexibility are desirable. First the ability to 
simulate a wide range of atmospheric pheno. cna is desirable. The ability to easily 
implement different microburst events from data sets such as JAWS is required. 

Another type of flexibility involves the freedom of the pilot to perform any 
maneuver desired. Many simulators assume a constant airspeed and commit the pilot 
to level or glide slope flight. The ability to go around in severe shears is denied. In 
microburst flights, the airspeed changes can be ±30 m/sec or more. The ability to go 
around or to test various escape strategies is highly desirable. 

The flexibility should extend to the ability to apply the wind and gust model 
to a wide variety of different aircraft. All frequencies of interest to aircraft response 
should be in the model. The ability to do this for aircraft of a wide variety of sizes and 
characteristics implies a nondimensional simulation. The concept of nondimensional 
simulation is explained in a later section. 

4. Easy implementation: Easy implementation implies code clarity, and 
simplicity, and portability. Portability means the ability to transfer code or data from 
one computer to another. If these attributes are missing, the method will not be 
accepted by the aviation community. 

While the above list is by no means exhaustive, it does include the major desirable 
attributes of a gust and wind model. 

In creating a simulation, the simulator is in effect creating a world of his own. 
This shadow world created for engineering purposes should contain enough spectral 
information so that all factors affecting the phenomenon under study are available. 

The following sections describe the creation of a turbulence and wind shear “world” 
which varies in space and time in a realistic manner. In the author’s opinion, the differ- 
ences in the realism of the present model and previously used one-dimensional models 
are similar to the differences in the creatures inhabiting our three-dimensional world 
and the one- and two-dimensional inhabitants of Lineland and Flatland [23]. 
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Creating realism for the sake of realism is not the function of engineering simulation, 
rather the engineering simulation is like a vignette. Relevant features are clear, crisp, 
and mathematically precise, while features not affecting the simulation are not. 

The model described in this document has many of the desirable attributes. Its 
greatest strengths are realism and flexibility. In these two areas, no currently available 
wind simulation model can match it. Easy implementation was a goal of this effort. 
Programs are in FORTRAN and data are stored in easily transportable integer formats. 
The one possible weakness of the model is in computational speed. The spatial model 
achieves maximum realism and flexibility but pays some price in speed. Nevertheless, 
for systems with enough central memory to store both the JAWS data and turbulence, 
speed should be sufficient to do real time simulations. Some minicomputers satisfy this 
requirement. 



CHAPTER V. OVERVIEW OF THE SPATIAL MODEL 


To this point, the aviation hazard posed by microburst-related wind shear was dis- 
cussed along with the characteristics of the JAWS data sets. Certain aspects of homogeneous 
turbulence theory relevant to Monte Carlo turbulence simulation were presented. Previously 
reported turbulence simulation methods and desired characteristics of wind model were dis- 
cussed. In this chapter all of these threads are pulled together to weave a fabric of computer 
created reality for the purposes o r aircraft design and pilot training and ultimately to save lives. 

Creating the final tapestry which achieves the ultimate reality as far as the wind 
environment is concerned is not a simple problem. The currently proposed technique is 
an attempt to combine the best available wind shear data set with three-dimensional 
simulated turbulence. The mixture of reality with simulated reality is the best that 
can be done within the current state-of-the-art. 

The need for the addition of turbulence is documented in this chapter. A recipe 
for combining measured data and JAWS data is explained, along with the method of 
implementation. The proposed technique is related to some similar one-dimensional 
approaches to data analysis and relevant concepts from these simpler cases are intro- 
du.ed. Finally, while the JAWS data sets are the best available in the world, they lack 
some infonnation which is necessary for realistic flight simulation. Methods for deriving 
the required parameters from aircraft measurements are presented. 

A typical grid spacing for the JAWS data is 200 tn whereas a desired grid spacing 
for the calculation of aerodynamic moments is 10 meters. In Figure 21 a comparison 
is made between the JAWS grid size and various aircraft. The box surrounding the 
planes represents vertical plane grid spacing for the July 14 case (200 m x 150 m). 

Figure 22 shows a desired grid spacing (10 m) compared to the Boeing 747 aircraft. 
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The spatial sampling frequency y s is 1/200 m'* and the resulting Nyquist spatial 
frequency is v^y = - 1/400 in'* . For an aircraft traveling at 80 m/sec this corre- 

sponds to a Nyqiist spatial frequency of fjqy = 0.2 Hz. Typically, the frequency of 
maximum aircraft response occurs at about 0. 1 Hz. According to Etkin [11], struc- 
tural response of an aircraft extends over a spatial frequency range from about 10"^ 
in'* to 10"“ m'*. Short period response occurs over a range from about 10'^ m'* to 10"^ 
in'*. Phugoid response cuts off at less than 10"^ m'*. The JAWS data contain Phugoid 
response frequencies, but only part of the short period response frequencies, and even less 
of the structural response frequency range. With regard to frequencies (airspeed = 80 m/sec), 
the JAWS cutoff frequency was 0.20 Hz while structural response goes up to near 10 Hz. 
Some modelers include structural bending modes in their simulators because they feel these 
modes add the correct “feel” to the turbulence. An aircraft structure encountering 
turbulence has a ringing response [24], The conclusion is that since structural and 
short term response of aircraft to turbulence are important for realistic simulation, and 
since the complete frequency range is not contained in the JAWS data, the high fre- 
quency turbulence must be added. 

The obvious question arising from the above discussion is, “how can turbulence 
be realistically added to JAWS and other data sets?” The addition cf these high fre- 
quencies are subject to several constraints. First, the model chosen is restricted to using 
mainly information available from the data set itself. For maximum realism, as little 
of the required information as possible should be generated from “rules of thumb” or 
from model equations. The available information includes the three velocity components 
(over the low frequency range) and spectral width which can be related to a gust 
standard deviation. Based on this available information the model to be discussed below 
was developed. The “glue” that holds the model together is the following equation. 


Uj(x,y,z,t) = Uj(x,y,z,t) + o^x.y.z.t) w^x.y.z) 


(38) 
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where Uj(x,y,z,t) arc the simulated winds, Uj(x,y,z,t) are the low frequency “smoothed” 
winds, aj(x,y,z,t) are the gust standard deviations, and Wj(x,y,z) are the zero mean, unit 
standard deviation, frozen turbulent velocities. 

The model defined in the above equation is written m its most general form but 
is specifically tailored for JAWS-type data sets. Uj are the three components of low 
frequency wind contained in the data sets, Oj, while written in vector form, in practice 
is isotropic and a scalar, o. It can be derived from JAWS second moment information. 
Notice that uj and Oj are written ->s functions of three spatial coordinates and time. 

As of this writing, time has not been included in the JAWS data sets. The reason for 
this is that the time variation may not be important for flight simulation because of the 
rapid aircraft transit times of wind shear phenomena. For this reason, the JAWS Project 
has concentrated its resources on other matters. 

The final term in Equation (38) is the frozen turbulence data base; frozen because 
Wj is not a function of time. A three-dimensional block of Monte-Carlo simulated tur- 
bulence is created and effectively stacked inside the JAWS data set in order to add the 
necessary small length scale phenomena to the coarser gridd^d JAWS data. In the current 
method the three-dimensional block is assumed to be isotropic turbulence. There is 
some basis for this assumption, since much of the anisotropy of the winds will be con- 
tained in Uj. Wj is assumed isotropic for convenience. The Monte Carlo simulation of 
nonisotropic turbulence is not far advanced in the one-dimensional case. In the three- 
dimensional case it is nearly nonexistant. Etkin, in his classic text [23], recommends for 
low altitudes (anisotropic turbulence) the three isotropic one-dimensional spectral func- 
tions be used with corresponding gust intensities. This is certainly a practical approach 
in the one-dimensional case, but the corresponding three-dimensional spectra-functions 
are not available. At least one group of investigators [25] have developed axisymmetric 
three-dimensional spectrum functions. These functions may prove useful in the near 
future, but for the present an isotropic model was used. 
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A question arising out of this discussion is, “Why use a three-dimensional turbu- 
lence model?” The reason is that by generating three-dimensional turbulence, all lateral, 
vertical, and longitudinal correlations are included in the data. Implicitly, when a one- 
dimensional simulation is done, the assumption is made that the aircraft is a point 
immersed in turbulence as is depicted in the left half of Figure 23. Turbulence is more 
complicated than this simple picture and two-dimensional turbulence is depicted in the 
right half of the figure. In fact, the spanwise variation of gusts can be quite high and 
gust differences measured in the B-57B Gust Gradient Program exceeded 10 m/sec (20 
kts). The B-57B has a 20 meter wing span and its size is compared to some wide-bodied 
transport aircraft in Figure 21. The Boeing 747 has a wing span of 60 meters. The 
impact of spanwise variation of gusts is undei study, but intuitively it seems that the 
variations have a significant impact on aircraft response. 



Figure 23. Assumptions of turbulence simulation. 

In Chapter II a cross-spectral model was developed. The block of three- 
dimensional turbulence contains the cross-spectral model as a natural subset of its 
properties assuming the three-dimensional von Karman spectrum is used to generate the 


turbulence. 
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Tlie model as proposed in Equation (38) was :n part inspired by a series of papers 
by Mark, and by Mark and Fischer [26-301. In these papers a one-dimensiona’ >-'rnu- 
lenee model defined by the following equation was investigated. 

w(t) = w s (t) + w t {t) = v s (t) + Oj<t) z(t) , (39) 

where w s (t) is the slowly varying part of an aircraft-measured velocity, Wjft) is the high 
frequency rapidly varying part of the velocity trace, z(t) is a Gaussian unit variance 
process (von Kartnan one-dimensional spectrum), and w s (t), a^t). and zG) are mutually 
independent processes. 

Mark points out that frequently a smoother knee than expected is measured in 
atmospheric turbulence spectra. He attributes this to the effect of a slowly varying gust 
intensity Of(t). The product op<t)z(t ) corresponds to a convolution in the frequency 
domain and if Of is a well-behaved function (no spikes in its spectrum), then the 
spectrum of the product will be a smoothed version of the spectrum of z. The result is 
a rounded knee in the spectrum of the product. 

Mark (1981) shows the •‘fleet on the autocorrelation function made up of parts 

as indicated in Equation (39), and presents data to support his model 129]. Figures 

24 and 25 are drawn from this report. In Figure 24, o w is the standard derivation of 

s 

the slowly varying part. wjt). The high frequency portion, w^t), contributes a rapidly 
decaying term to the autocorrelation and w s (t) contributes a slowly decaying term. In 
Mark’s model, w s (t), Ojtt), and z(t) are presumed independent processes. In Figure 25 
a corresponding measurement is presented. In the frequency domain the role is reversed. 
The spectrum of w s (t) will decay rapidly while the spectrum of Wjtt) decays slowly. 

The summary of the analyses by Mark and by Mark and Fischer was for or.e- 
dimensional turbulence measurements, but most of the results carry over to the three- 
dimensional model of Equation (38). A notable exception concerns the mutual inde- 
pendence of uj, Oj, and Wj. The mutual independence assumption was made by Mark 
to facilitate his analysis, but in fact u^ and Oj should be related. We expect this because 
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Figure 25. Auto-correlation function of vertical veiocity record in 
mountain wave conditions (airspeed 197 m/sec). 

of mechanical production terms in Reynold’s equations involving gradients in Uj. The 
relationships may not show up in linear statistical relationships, however. In any case, 
the independence assumption is not required in the present model. 

The terms in Equation (38) were discussed above, and now a discussion concern- 
ing the assembly of the three terms into a wind model is presented. For the JAWS 
July 14 data set the grid consisted of 60 x 60 x 11 points. At each of these points is a 
measured value of the three velocity components, and spectral width. This corresponds 
to 158,400 data points. The' points were measured on a 200 m x 200 in x 150 m grid. 
A desirable spectral resolution for the turbulence is 10 m x 10 m x 10 m. Generation of 
a block of turbulence 12 km x 12 km x 1.5 km with 10 m resolution corresponds to 
6.5 x 10^ data points and is not feasible. Instead, a relatively small block of turbulence 
can be generated and effectively stacked to fill the JAWS volume. Stacking is equivalent 
to moving the block around within the JAWS volume as the plane flies through and 
begins to leave the turbulence. The situation brings to mind a small boy playing with 
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a caterpillar. As the caterpillar hangs over the edge of his hand the boy obligingly offers 
his other hand and the worm’s journey of exploration continues. In his travel across 
each hand the caterpillar is unlikely to retrace the exact path of his previous journey 
because he begins his journey at a slightly different point each time and travels with a 
slightly different heading. The cycle continues as long as the boy desires and the result 
is that the caterpillar has his exercise and the boy has the pleasure of his company. 

The periodic shifting of the data base is equivalent to stacking blocks of turbu- 
lence within the block of JAWS data. Two ways of stacking the blocks come to mind. 

As the aircraft of the simulated flight passes out one side of the block its motion can be 
reflected back into the block of turbulence This procedure is similar to video games in 
which the electronic ball bounces off the wall (angle of incidence equals angle of reflec- 
tion). This method, through complicated has the advantage that the aircraft encounters 
no abrupt discontinuities in the turbulence field. This approach is equivalent to stacking 
two types of blocks in a special pattern (see Figure 26). The two types of blocks are 
reflections of each other, the same as left and right hands. The blocks are stacked so 
that similar sides touch each other. 

Computationally a simpler approach is to stack one type of block on itself with 
the same orientation. Position within the block is calculated using congruence arithmetic. 

X T = X mod (X Tmax ) 

Y T = Y mod (Y Tmax ) (40) 

Z T = Z mod (Zxmax) • 

In this equation, Xy is the X-position within the block of turbulence, X is the real X- 
position in space, Xy max is the maximum extent of the frozen turbulence, and so on. 
This expression means, divide X by Xy max , determine the remainder which is Xy. If 
X = 6500 m and Xy max = 600 m then Xy = 500 m. The value of Xy always lies 
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Figure 26. Stacking blocks of frozen turbulence. 
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between 0 and Xy max . This method is equivalent to video games in which the electronic 
“projectile” disappear from one side of the screen and reappears instantaneously on the 
opposite side. An anticipated disadvantage of this technique is a discontinuity in 
the turbulence as the plane passes from one side of the block to the other. In practice, 
the discontinuity is small or nonexistent so the computationally simpler technique of 
Equation (40) can be used. 

One caution on the use of the above equation should be discussed at this time. 
The frozen turbulence data base generated for use with JAWS data is nondimensional 
in space. Transformation of a complete block of turbulence must be achieved at one 
time with one length scale of turbulence. At the same time, the turbulence length scale 
should vary through the simulated atmosphere. The most commonly studied variation is 
with height above ground. The presence of microbursts in the planetary boundary layer 
(PBL) also creates a lateral variation of length scale in the atmosphere. Convective 
storms imbedded in planetary scale flows have been observed to create an obstruction to 
the larger scale flow and even shed vortices. Apparently convective obstructions affect 
turbulent scale lengths in the atmosphere. There is every reason to expect that micro- 
bursts which are embedded in larger scale How in the PBL like chocolate candy kisses 
will similarly affect length scales in the PBL. In order to account for desired changes in 
this length scale, the turbulence is generated in dimensionless space. In dimensionless 
space the block size is constant, but in dimensional space the block size changes 
isotropically with the length scale. In the space domain an increasing turbulent length 
scale corresponds to increasing block size, and in the spatial frequency domain a decreas- 
ing block size. The problem with indiscriminant use of Equation (40) is that the blocks 
are stacked on top of each other and as each contracts or expands the whole stack 
similarly contracts or expands. If the stack contracts along its lower left hand comer 
the net motion at the right side of the stack is multiplied by the number of blocks in 
between. The result is an unreal change in aircraft position relative to the turbulence. 
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Care must be taken that stack contraction is about the location of the aircraft and this 
is easily implemented. The nondimensional generation of turbulence mentioned above 
will be described in detail in the following chapter. 

The preceding discussion brings us to the one parameter not available from the 
JAWS data set, i.e., the length scale of turbulence. The distribution of this parameter 
must be either modeled, or derived from aircraft experiments or from laboratory 
measurements. This desired length scale is not the overall length scale but rather the 
length scale associated with wj(x,y,z) in Equation (38). 

Interestingly, some evidence suggests that the spanwise length scale increases 
through a microburst {31]. The relevant information is presented in Figure 27. Shown 
are the three velocity components measured at the center boom of the B-57B aircraft 
during a suspected microburst encounter. The encounter occurred about 78 seconds 
into the run when a sudden sharp headwind increase was observed. This increase 
corresponds to the sudden decrease in the longitudinal velocity component, u. During 
this period the B-57B encountered a 15 m/sec headwind increase over a distance of about 
130 m. This drastic windspeed change is believed to occur as the aircraft crossed a 
microburst front. Aircraft altitude was about 400 m AGL and from Chapter II we recall 
that the typical outflow depth is 600 m. The increasing headwind was followed by a 
decreasing headwind and then a tailwind all over a period of about 20 seconds. The 
feature had an extent on the order of 2 km which is consistent with the microburst 
hypothesis. The horizontal wind vector shifts were associated with a strong (10 m/sec) 
downdraft which reinforces the microburst idea. 

The interesting part of the figure is the lower three graphs which depict wingtip 
to wingtip velocity differences. The amplitudes of the velocity difference traces increase 
gradually from the start of the run but then decrease in all three components through 
the suspected microburst. This unexpected result can arise from three possible sources: 
(1) the gust intensity decreases through the microburst, (2) the lateral length scales 
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increase through the microburst, or (3) a combination of (1) and (2). From the traces 
at the top of the figure, it is difficult to believe that an intensity decrease is the sole 
source of the velocity difference decrease. The conclusion then would be that a length 
scale increase occurs through the microburst. Some caution must be exercised in drawing 
conclusions from one flight through a microburst, but this type of information is vital 
for modeling distributions of L. 

The problem with obtaining L distribution from aircraft measurements is the 
difficulty of finding microbursts to fly through even with Doppler radars to guide the 
aircraft. At the same time, more information is required. Laboratory experiments may 
help provide the much needed information. Means of obtaining this information in the 
laboratory are discussed in Appendix A. 

The dearth of knowledge on L distributions is similar to the situation turbulence 
modelers found themselves in prior to the development of two-equation turbulence 
models [32]. Mixing length distributions had to be specified based on experience. 

Finally two-equation models were developed which included a partial differential equa- 
wion which effectively specified distributions of the mixing length. The close analogy 
with the present problem suggests that a similar approach to the present problem might 
be attempted. The most widely utilized two-equation model is the k-e model which 

7 

contains transport equations for turbulent kinetic energy k = 1/2 (u~ + + w~) 

and turbulent kinetic energy dissipation e. The equations for this model are presented 
by Launder and Spalding [33] in the following form. 


Dk 

Dt 
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p dx^ 


^eff 3k \ 

^°k 3x k J 




(41) 
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Dt 



C 1 **eff e 
p k 




where M e ff = Qipk^/e and recommended values for curve fit “constants” are 
= 0.09, Cj = 1.44, C 2 = 1.92, a k = 1.0, o c = 1.3 . 

The length scale of turbulence is defined by 

L - C D k 3 % (43) 

Here Cq is another constant of proportionality. 

Observe that if JAWS estimates of k = 3/2 o ^ based on spectral widths are used 
then Equation (40) need not be solved. If the Uj in Equations (41) nd (42) are assumed to 
be the JAWS velocities, appropriate difference forms of Equation (42) can be developed 
and with appropriate boundary conditions for e and with the aid of Equation (43), a distri- 
bution of L could be calculated. In reality the “constants” in the above equations vary from 
one flow to another and would need to be defined for microburst type flows. The above 
equations were developed by making certain assumptions which may not be appropriate 
for the JAWS microburst situations. Nevertheless, the above provides an engineering 
approach to a difficult problem. The method can be refined and “tuned” if experimental 
data are available to tie down the parameters in the above equations. 

Flight simulators would use the spatial model as follows. The frozen turbulence 
data base would be giver, to the particular simulation group on magnetic tape. The data 
on tape would be transferred to high speed, random access mass storage, e.g., disc storage. 
The size of the data base provided would depend upon the available mass storage for the 
particular computer. For reasons given in the next section, it is desirable to generate the 
data base in one large block. The total block may be too large for many simulators, and it 
is not necessary to use the entire data base. The main block can be divided into subblocks 
which will fit into available mass storage. If additional mass storage became available at a 
later time, then contiguous blocks of data could be shipped to create a larger data base at 


the installation. 
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Finally, Figure 28 shows how drastically the JAWS data truncates assumed 
spectra. This last figure in this chapter serves to reemphasize some points made herein. 
First, for realistic wind modeling for flight simulation, turbulence must be added to the 
JAWS data. Secondly, the model described by Equation (38) is a model based on 
necessity. It is not a scientifically accurate turbulence model, but it is a realistic 
engineering approach to a difficult problem. It uses all available information from JAWS 
and requires only turbulence length scale distributions for completeness. While the 
frozen turbulence is Gaussian, the winds generated by the model of Equation (38) are 
not, as real atmospheric winds are not. The model contains all three-dimensional 
correlations and spectra as any realistic model should. Finally, the model was specifically 
tailored for use with the JAWS data base, but is not restricted to it. If Oj and Uj were 
known for winds on Jupiter, it is anticipated that flight vehicle entry into the Jovian 
atmosphere could be simulated with the proposed model. 
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Figure 28. Cutoff frequencies for the JAWS data sets (von Karnian spectra). 
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CHAPTER VI. FROZEN TURBULENCE GENERATION 
A. The Fast Fourier Transform Approach 


Three-dimensional turbulence generation was achieved using Fast Fourier Trans- 
form (FFT) techniques. Figure 29 is similar to Figure 20 but differs in that the filter 
function, input, and output are functions of three independent variables (four if time is 
included). If the FFT is used, the Dryden model has no computational advantage over 
the von Karman model. Hence, the von Karman model was used. Generation of turbu- 
lence using the FFT requires some awareness of certain FFT properties. These properties 
are discussed with some derivations given in the appendices. In this chapter all frequen- 
cies, spectra, and filter functions are assumed dimensionless. 


n i ( f 1< r 2< r 3) 


WIDE BAND 
GAUSSIAN NOISE 


H (v\, v 2, 1^3) 


Wj ( r 1- r 2 , r 3 ) 

HOMOGENEOUS, ISOTROPIC 

THREE-DIMENSIONAL 

TURBULENCE 


Figure 29. Generation of three-dimensional turbulence. 


Since the turbulence is generated digitally, sampling frequencies in the space and 
spatial frequency domain play an important role. The relationships are given by the 
following equations 
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where v s . is the sampling frequency corresponding to the rj direction, Mj is the number 
of grid points in the rj direction, Ar>j is the grid spacing in the spatial frequency domain 
corresponding to the rj direction, and Arj is the rj grid spacing in the space donum. 

The FFT is a fast computational implementation of the Discrete Fourier Trans- 
form (DFT). One-, two-, and three-dimensional DFTs and their inverses are defined in 
the following equations. 

M-l 

X k = £ x n e'J 2,rnk / M k =0, 1 M-l (46) 

n=0 


M-l 

x =- V X k e> 2 ™ k / M n = 0, 1 M-l (47) 

11 M 

k=0 


X kj,k 2 


Mj-1 M 2 -l 

^-j2rr(njkj/N] + n^k^/N^) 

n j =0 n 2 =0 



kj = 0, 1 Mj-1 

(48) 


M r l M 2 -l 


1 > n 2 ~ M, M-i Xk l’ k 2 C 

1 *• k j =0 k 2 =0 


+j2ir(n jkj/M j + n->k-*/M 2 ) 


nj = 0, 1, . . . , Mj-1 
(49) 


M , — 1 M 2 -l M 3 -l 

X k.,k 2 ,k 3 = X S £ x nj,n 2 ,n 3 
nj=0 no=0 n 3 =0 

Mj-1 M 2 -l M 3 -l 
1 V V V V +j2n(n 1 k 1 /M ] +n,k,/M, + n,k 3 /M 3 ) 


-j2rr(njkj/Mj + n-)k->/M 2 + n 3 k 3 /M 3 ) 

kj = 0 1, ... , Mj-1 
(50) 


nj = 0, 1, . . ., Mj-1 
(51) 
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In Equations (46) through (51) n and lower case x refer to the space domain and k and 
upper case X to the frequency domain. In the above six equations, even numbered equa- 
tions, i.e., (46), (48), and (50) are the forward transforms and the odd numbered ones 
are the inverse transforms. The three different transform pairs are represented becaus -• in 
the discussion to follow symmetries associated with the one- and two-dimensional trans- 
forms show up in the three-dimensional DFT. These symmetries are of considerable 
importance to the generation of three-dimensional turbulence. 

The three-dimensional continuous ourier transform and its inverse are given by 
the following equations. Since all DFT implementations in this report are by means of 
the FFT, these two terms will be used interchangably 

X^i,!^) = J J j x(rj,r 2 ,r 3 ) e + X - V - + drj dr „ dr3 (52) 

—DO -OO -OO 

oo oo oc 

x(rj,r2.r 3 )= f f f X^j,^,^) + V - T ~ + dnj dv 2 di> 3 (53) 

-OO -C 5 -OO 

By comparison of (52) and (53) with (50) and (5 1 ), approximations for the Fourier 
transfoim and inverse are given by 

X(kjA^|, k 2 A^ 2 , k 3 Ar> 3 ) ~ Aq Ar 2 Ar 3 DFT (xfnjArj, 02 ^. t^Arj)] (54) 

x(njArj, n 3 Ar 3 ) ~ MjA^j NtaAi^ M 3 Ai> 3 DFT"^ [X(kjAvj, k 2 Ai> 2 , k 3 Ar> 3 )] 

(55) 

Similar expression^ apply for the one- and two-dimensional transforms. 

Figure 29 seems to imply that two three-dimensional transforms are necessary for 
turbulence generation. First, noise generated in the space domain is transformed to the 
frequency domain, multiplied by a filter and then the inverse transformed back to the space 
domain. Actually, only one transform is required because the noise can be generated in 
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the frequency domain since transformed noise statistical properties can be calculated. 
Important to this discussion are the symmetry properties of transformed real noise for 
the one, two-, and three-dimensional cases. Derivations follow easily from the transform 
definitions. These properties are presented in Table 2. 


1 . 

2 . 

3. 

4. 

5. 

6 . 

7. 

8 . 
9. 

10 . 

11 . 


Table 2. Symmetry Properties of Transformed One-, Two-, and 
Three-dimensional Digital Functions 


One-Dimensional Symmetry 


Expression 
X(M-k) = X*(k) 


Observation 
Im[X(M/2)] = 0 


Two-Dimensional Symmetry 

X(Mj-kj,0) = X*(kj,0) I m [X(Mj/2,0) = 0] 

X(0,M 2 -k 2 ) = X*(0,k 2 ) I m [X(0,M 2 /2) = 0] 

X(N r kj,N 2 -k 2 ) = X*(kj,k 2 ) I m [X(M j /2, M 2 /2) = 0) 


Three-Dimensional Symmetry 


X(Mj-kj,0,0) = X*(kj,0,0) 

X(0, M r * , 0) = X*(0,k 2 ,0) 
X(0,0,M,-k 3 ) = X*(0,0,k 3 ) 
X(M r k l5 M 2 -k 2 ,0) = X*(kj,k 2 ,0) 
X(M r kj, 0, M 3 -k 3 ) = X*(k 1 ,0,k 3 ) 

X(0, M 2 -k 2 , M 3 -k 3 ) = X*(0,k 2 ,k 3 ) 
XfMpkj, M 2 -k 2 , M 3 -k 3 ) = X*(kj,k 2 ,k 3 ) 


I m [X(M 1 /2,0,0)] =0 
l m [X(0,M 2 /2,0)l = 0 
I m [X(0,0,M 3 /2)l =0 
I m [X(M j/2, M 2 /2, 0)] =0 
I m [X(M,/2,0,M 3 /2)) =0 
I m [X(0, M 2 /2, M 3 /2)J = 0 
I m [X(Mj/2, M 2 /2, M 3 /2)J =0 
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The symmetry relations expressed in Table 2 seem quite complex, but all can 
be incorporated into a simple geometric interpretation. The geometric interpretation is 
reflection about the center point. In the one-dimensional case the reflection is about the 
center point (k = M/2) of the line, in the two-dimensional case, points along the two axes 
are reflected about their respective center points (kj = Mj/2, k-> = 0 ; kj= 0 , k-> = M 2 / 2 ) 
and in the grid interior about the center point of the interior (kj,k2 =£ 0, reflection 
about kj = Mj/2, k-> = M->/2). In the three-dimensional case one-, two-, and three- 
dimensional reflections must be made about the three axes, the three planes (kj-k2, 
kj-k3> k2-k3>, and about the interior symmetry point (kj.k2.k3) = (Mj/2,M2/2,M3/2) 
in the interior of frequency space. The Hermitian symmetries are depicted for the one- 
and two-dimensonal cases by Figures 30 and 31 . Showing the symmetries for the three- 
dimensional case is very difficult but the principle is the same. Figure 32 indicates the 
grid points at which the transform of a real function must take real values. If these 
points fall on an axis, Hermitian symmetry exists about the point on the axis. If the real 
point falls on the kj^, kj-k3, or k-,-k3 plane, Hermitian symmetry about the point 
exists in the plane, and similarly for the interior point. Each point other than these 
eight real points possesses Hermitian symmetry with another point in the grid through 
the real point and extended a distance equal to the distance from the original point to 



EXPRESSION 1, TABLE 2 



Figure 30 . Correspondence of complex conjugate pairs in the 
one-dimensional transform domain. 
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Figure 32. Points in the transform domain where the 
transforms must be real. 

the real point. The end of the line segment falls at the corresponding grid point. Values 
of the transform at the original point and the new point will be complex conjugates of 
each other. This fact is a consequence of the requirement for real turbulence, i.e., not 
complex. In one dimension, there are two real points (k = 0, M/2), in two dimensions, 
four real points, three dimensions, eight real points, and so on. 

Linearity is a well known property of the DFT and DFT^. Since the wide band 
random noise to be transformed is Gaussian, then the transformed noise is al. Gaussian. 
This being the case, if the expected value and variance of the transformed noi— were 
known, then all the information required for generating noise in the frequency domain 
would be known. The mean and variance of transformed noise for the one- and tliree- 
dimensional cases are calculated in Appendix D and summarized in Table 3. With the 
information available in Table 3, generation of three-dimensional transformed noise in 
the frequency domain is possible with a simple Gaussian random number generator. 
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Once generated, the noise must be rearranged according to the Hermitian symmetries of 
Table 2 so that when transformed back to the space domain the resulting turbulence 
takes real rather than complex values. 

Table 3. Expected Value and Variance of Transformed One-, and 
three-dimensional noise with zero mean value 


One-Dimensional Case 
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B. Derivation of Filter Functions 


A simple approach for deriving the one-dimensional filter function is presented 
and then generalized to the three-dimensional case. Since the turbulence is generated 
digitally, the frequency cutoff, i.e., the Nyquist frequency, is the sampling frequency 
divided by two O s /2). The DFT is an approximation to the continuous Fourier trans- 
form evaluated over these limits (see Figure 33). The magnitude of the constant 
spectrum was selected as o n 2 /i> s so that the following equation was satisfied. 


V 2 

/ *n dl, s = °n 2 

-V 2 


(56) 


The expression for the output spectrum for a linear filter in terms of the input 
spectrum is 

<h w 0') = \m H*(tO 4> n (»0 = |H(n)| 2 *„(»') (57) 

The output spectrum is the desired spectrum and is known while the input spectrum is 
as shown in Figure 33. The solution for H(n) is 


y/ v s*v, 
H(?) " 


ja(t') 

e 


(58) 


where ot(v) is any arbitrary phase function. For this study, a(v) = 0 was selected. 

The one-dimensional case can be generalized to any number of dimensions. For 
the three-dimensional case the analog to Equation (56) is 



~ v s\/2 


v d 2 v s3l 2 

f f *n 

'"S2 I 2 ~ v s3 / 2 


dv 2 dn 3 = o n 2 


(59) 
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This equation implies a value of 4> n of o n “/(^ s i ^2 ^ 3 ). Equation (57) still applies 
with each term being a function of these variables. The new filter function is given by 

^s2 ^3 


n^i PiPj) = — 





\ 






-f % 


Figure 33. Input noise spectrum. 

This equation is the general form of the Filter function. $ w for this study is the 
three-dimensional von Karman spectral function for either the u, v, or w velocity com- 
ponent. In Chapter II the dimensional forms of the von Karman spectra were presented. 
In Chapter V the need for nondimensional spectra was described. The dimensional 
spectra are transformed by = £jL where is the dimensional wave number. Then the 

A 

dimensionless spectrum 4> w is given in terms of the dimensional spectrum 4> w by the 
following relation: 

a (*\ v 2 "3 

♦w^i’W) Av i Av 2 dv 3 = ^wV l ’ r ’ r 

With this transformation, the nondimensional spectra given in Equation (13) become 

*ii _ 440 tt 3 a4 - v ? > 

^ [1 + (27rai>)^l ^/6 


) 


d^j d^ d^ 


(61) 


( 62 ) 
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With Equations (60) and (62) and the results of Tables 2 and 3 all the tools are available 
for generation of noise in the frequency domain. The sequence is as follows. Generate 
three-dimensional noise over half of three-dimensional space then fill the rest of the space 
using the Hermitian symmetry properties of Table 2. Make sure the noise variances 
satisfy the restrictions of Table 3. Multiply the transformed noise point by point by a 
sampled version of the dimensionless filter function derived from Equations (60) and (62) 
above. Apply the inverse FFT and the result in the space domain is three-dimensional 
dimensionless turbulence which can be added to mean wind data sets. Codes for per- 
forming each of these steps are described and listed in the appendices. 

The generated turbulence contains all the correlations with lags in each of the 
three space directions. This being the case, if a line of turbulence, say in the rj direction 
is selected and the one-dimensional spectrum calculated with an FFT the result should be 
an approximation to the one-dimensional continuous spectrum. To check tills, consider 
the generated turbulence given by 


w i 


M r l 

M-,-1 

M 3-1 




An 3 £ 

t 

E 

H(kj,k 

2> k 3> 

N(k!,k 2 ,k 3 ) 

kpO 

k 2 -0 

o 

II 







( 

n 1 k i 

n 2 k 2 n 3 k 



exp 

+j27rl 

Mj 

+ ■ + 

m 2 m 3 


Then wfn^n^.n^) is transformed with the following equation. 


(63) 
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V" -j2irni£i/Mi 

w(lj,n 2 ,n 3 ) = Arj > wOij.no^) e 1 1 1 
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(64) 


If Equation (63) is substituted into (64), and algebra performed as in Appendix E, the 
final result is 
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m 2 -i m 3 -i 

® w(ll)= zk^ £ ^ > n^ 1 l ,k 2 ,k 3^ ^ v 2 ^ v 3 (65) 

1 k 2 =0 k 3 =0 

This is the digital analog to Equation (17). The factor l/A^j must be present because 
© w is an energy density, i.e., the turbulence energy between frequencies Vj and + dvj 
is given by 0 w (r>j) d^j. Equation (65) is consistent with previous derivations. 

A description of the means for turbulence generation was presented in this 
chapter. The technique is a three step procedure, namely (1) generation of three- 
dimensional noise, (2) multiplication by the sampled spectral function, and (3) inverse 
transformation by FFT to the space domain. An attempt was made to help the reader 
understand the important aspects of the problem by presenting geometrical interpreta- 
tions of the procedure. A clear understanding of these interpretations helps in under- 
standing the codes used for the turbulence generation. Two FORTRAN programs were 
written, one for the first two steps above, and another for the final step. These pro- 
grams are described in Appendices F and G, respectively. 


- 





CHAPTER VII. USE OF WINDSHEAR DATA SETS WITH 
SIMULATED TURBULENCE 


A description of the addition of turbulence to a JAWS data set is presented. The 
JAWS July 14 case was selected because some aircraft data from the same day were 
available for estimating gust intensities and length scales of turbulence. While these 
aircraft data were not from the Doppler radar measurement region, the results should be 
applicable in a general way to the aircraft measurements. This conclusion was based on 
the fact that July 14 was a day when many microbursts occurred over a large area. 

Any gross variation in turbulence characteristics from one micioburst group to another 
on this day seems unlikely. 

The JAWS data set characteristics were described in Chapter II and Table 1 
indicates that grid spacing on the July 14 case was 200 m by 200 m by 150 m. For 
simulated flights using this data, interpolation of the JAWS winds was required. The 
interpolation procedure was based on bilinear Lagrange polynomial basis functions. A 
good description of two-dimensional interpolation is presented by Prenter, 1975 [34]. 

The two-dimensional method is easily generalized to three or n dimensions. For three 
dimensions, the general interpolation form is 

M , — 1 M 2 -l M 3 -l 

U(x,y,z) - ^ ^ ^ ^ni,n 7 ,n 3 ^ni,n 7 ,n 3 ^ x, y ,z ^ > (66) 

nj=0 n 2 =0 n 3 =0 

where U(x,y,z) is the interpolated quantity, U n ^ n is the value of the quantity at 
the grid points and is a group of basis functions, one for each grid 
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point. The simplest set of basis functions for our need is the set of bilinear Lagrange 
polynomials. Each of these basis functions is continuous so any finite sum of the 
functions is also continuous. Therefore, continuity of the interpolated quantity U(x,y,z) 
is assured. 

The Lagrange basis functions P n are defined at each grid point so that 

P n 1 ,n 9 ,n 3 (x n 1 ’ym’ z n 3 ) = Xn^nv and z n 3 are x - and z coordinates of grid 

point (nj,n 2 ,n 3 ). At all surrounding grid points P n m (x m ’ y nv z n3 = °- ln 

1’ 2’ j 1 L 3 

practice the summation of Equation (66) is only over the surrounding eight grid points. 
In equation form the basis functions associated with a particular cell are defined by 


Pj(x c ,y c ,z c ) = (Ax-x c ) (Ay-y c ) (Az-z c )/(Ax Ay Az) (67) 

P 2 (x c ,y c ,z c ) = x c (Ay-y c ) (Az-z c )/(Ax Ay Az) • (68) 

P3( x c>y c > z c) = x c yc (Az-z c )/(Ax Ay Az) (69) 

P4(x c ,y c ,z c ) = (Ax-x c ) y c (Az-z c )/(Ax Ay Az) (70) 

P5( x c >y C ’ z c) = ( Ax_x c ) < A y-y c ) z c ^ Ax A y Az ) (7i) 

p 6^ x cyc’ z c) = x c ( A y-y c ) z J (Ax A y Az ) (72 ) 

P 7 (x c ,y c) z c ) = x c y c z c /(Ax Ay Az) (73) 

p 8( x c>yc’ z c> = ( Ax-x c ) yc z c/ (Ax A y Az > > (74) 


where x c ,y c ,z c are the coordinates within the cell as defined in Figure 34; Ax,Ay,Az are 
the x, y, and z grid spacing, respectively; and Pj is the basis function corresponding to 
grid point i as defined in Figure 34. 

Notice that gridpoint one lies at (x c ,y c ,z c ) = (0,0,0), grid point two at (x c ,y c ,z c ) = 
(Ax, 0,0), etc. Notice also that Pj is equal to 1 at grid point i, is always nonnegative, and 
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Figure 34. Definition of interpolation variables. 

is zero at surrounding grid points. Outside of the cell shown in Figure 34, each of the 
eight basis functions have different definitions, but these alternate definitions are of no 
concern to the interpolation of U within the cell. 

Interpolation within the grid of Figure 34 is achieved with the following equation. 

8 

U( x,y,z) = £ Uj Pj(x,y,z) , (75) 

i-1 

where Uj is the value of 0 as the ith grid point. 

The material in this chapter was implemented for the most part by the computer 
program listed in Appendix G. This program was written on a minicomputer with 
limited central memory. As a result, data was swapped from central memory to disc 
files and back. Three-dimensional arrays could not be conveniently stored on disc. The 
arrays had to be stored in linear records. For this reason, the JAWS data was stored 
cell by cell on disc. In other words, the eight values of U at each of the eight grid 
points of Figure 34 were stored contiguously. An adjacent cell would have four of the 
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same values stored so that memory-wise, this procedure is inefficient. Computation- 
wise it is efficient because accessing the grid values of the cell requires only one READ 
statement as opposed to four required by another procedure. More of these details arc 
described in Appendix G. 

The addition of turbulence to the JAWS data required information on the varia- 
tion of gust intensity and turbulent length scale. At the time of writing, JAWS second 
moment information was not available. An equation for the distribution of gust intensity 
was necessary for simulation. As a guide, consider Figure 27. This figure shows veloci- 
ties measured by NASA’s B-57B aircraft during flight through an apparent microburst 
on the day after the July 14 JAWS case. An interpretation of Figure 27 was presented 
in Chapter V. Evidence of decreased gust intensity and increased turbulent length scale 
were described. For the purposes of demonstration a combination of both was assumed. 

An exponential type decrease through the JAWS microburst was assumed. Lateral 
variation was written in terms of distance r from the microburst. Moving away from the 
microburst center, gust intensity increased to a maximum and then decreased to zero 
exponentially at great distances. A vertical exponential factor was added to make 
intensity a maximum at 300 m AGL. The functional form is given by 

ff = K°max “ (°max/2) exp"t (r /°- 75 ^ ' ) exp(-r/50) exp ((z-300/200)]. (76) 

where r is horizontal distance from the microburst center in km and z is the altitude in 
meters. 

A simple vertical relation for turbulent length scale was selected from Reference 
35. The functional form is 

L = 31.5 (z/18.3) 0 c4 meters (77) 

Turbulence generated for use with the JAWS data consisted of a 64 x 64 x 64 
array generated at dimensionless frequencies of = 50, ^ - 50, = 50. One concern 
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was the accumulation of error in the calculation of an FFT this large. The three- 
dimensional FFT does not have an error problem, at least for arrays of this size. Test 
cases were done and error was on the order of 10'^ which is the same as the word 
precision (32 bits per word, 24 bits for the mantissa). 

Interpolation of the turbulence between grid points was required. Tatom and 
Smith, '982 [36] offer some guidance on this point. They show that zeroth order 
interpolation (stair step turbulence) causes aliasing to exactly cancel the effect of digi- 
tization of the data as long as the turbulence is sampled at frequencies the same as or 
lower than generation frequencies. This procedure was followed. 

To illustrate model performances a point “airplane” was flown at constant 
ground speed through the data. The program in Appendix G that implements this model 
has a pair of “knobs”, one of which turns the turbulence up or down, and the other 
turns the wind shear on or off. Hence, the JAWS data without turbulence, the turbu- 
lence generated without JAWS data, or JAWS data with turbulence can be examined. 

The point “airplane” path was selected J o begin at a point so that a true heading 
of 35° carried the plane through the niicroburst center. A simulated ILS approach w.th 
three degree glide slope was flown beginning at an altitude of 350 meters. 

Figure 35 shows the east-west velocity component encountered in flight without 
added turbulence. Figure 36 is the same flight path with turbulence added. Notice 
that the length scales are decreasing toward the end of the flight as the plane approaches 
the ground. The turbulence in this figure appeals quite realistic (compare with Figure 
27). 

The traces of turbulence only are plotted in Figure 37. The data was sampled 
at a very high rate so that each data point was shown. Input parameters were arranged 
for horizontal flight through the block in the east-west direction. The middle curve is a 
plot of the points ix = 1-64, iy = 1, iz = 1, the upper curve ix = 1-64, iy = 2, iz = I, and 
the lower curve, ix = 1-64, iy = 32, iz = 32. The first two adjacent curves are highly 
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Figure 35. Lasi-west JAWS velocity on a simulated ILS 
approach (three degree glide slope). 
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Figure 36. East-west JAWS velocity component plus turbulence on 
a simulated ILS approach (tliree degree glide slope). 
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Figure 37. Comparison of lines of turbulence from the simulated turbulence. 

correlated. The third curve shows only slight correlation with the first two, as was 
expected. These three traces show the kind of transverse (lateral and vertical) three- 
dimensional correlation which was missing from previous one-dimensional turbulence 
simulations. 

The storage requirements for this realistic three-dimensional simulation are not 
as large as one might think. The July 14 JAWS case was provided on a 61 x 61 x 11 
grid. For three velocity components and one second moment, the number of words is 
163,724. The turbulence was stored on a 64 x 64 x 64 grid and with three components 
the number of words is 262,144. The JAWS and the turbulence data were stored in 
integer form in only two bytes. The total storage requirement is only slightly greater 
than 850 kilobytes. The Perkin-Elmer 3250 minicomputer can be outfitted with 32 
megabytes of central core storage most of which is available to the user. The conclu- 
sion is that very sophisticated simulations can be performed on small computers. 



CHAPTER VIII. SUMM ..Y AND CONCLUSION 


The characteristics of the spatial model can be summarized as follows: 

1. The spatial model contains three components of real wind shear varying over 
space with corresponding exponents of simrlated turbulence also varying over the 
three space dimensions. Bee ise wind and gust variation over the body of an aircraft 
are available, all aerodynamic loads and moments can be calculated. 

2. The simulated turbulence is nonlinear, non-Gaussian, and conforms to the 
von Karman three-dimensional spectral model. 

3. Because of the conformance to the three-dimensional spectr-, the model 
contains cross-spectral information for each component. 

4. By virtue of its three dimensionality, the spatial model permits flight simula- 
tions of any maneuver without diminishing the validity of the simulation. 

5. The model is highly flexible. Any flow field about which ensemble average 
velocities, gust intensities, and turbulent length scales are known can be simulated. 

The spatial model was implemented on a Hewlett-Packard F series minicomputer. 
The three-dimensional turbulence was generated in a 64 x 64 x 64 block and showed the 
three-dimensional correlation expected of it. The resulting turbulence when added to 
the JAWS data resulted in a simulated wind trace virtually indistinguishable to the eye 
from real data measured with the B-57B aircraft. 

An extension of the Houbolt-Sen cross-spectral model was presented and com- 
pared with data from the Gust Gradient Program. The theoretical curves showed some 
variation from the measured data. This variation is believed to be a probe effect, aliasing, an 
effect of the nonfrozen nature of the turbulence, or a combination of the three factors. 
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The spatial model requires one statistic not available from the JAWS data; the length 
scale of turbulence, L. The effect of microbursts on distributions of L should be a subject 
for future research. This research should include flight and laboratory measurements. 

Flight measurements are required because of the inability to model atmospheric scale 
Reynolds number in the laboratory. Laboratory measurements are required because of the 
difficulty of locating and flying through a microburst. Although laboratory measurements 
cannot provide quantitative results, they should provide trends. 

Laboratory measurements can contribute significantly to the understanding of the 
effect of topography on microbursts. The primary controlling parameter for microburst 
shape is the Froude number which can be modeled very well in the laboratory. Topographic 
effects in the JAWS data should be identified and carefully considered where the data are 
used with the spatial model. 
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APPENDIX A. MICROBURST RESEARCH IN THE LABORATORY 


The dimensionless parameters governing microburst-type flows are the Reynolds and 
Froude numbers (Re = VD fv, Fr = V\/g(j3AT)D), where 0 is the volumetric expansion 
coefficient. For real-atmosphere microbursts typical velocity and length scales are 10 m/sec 
(downdraft) and 600 m (depth of out-flow). The atmospheric Reynolds number is then on 
the order of 4 x i0^. A typical microburst AT is 3 deg C which implies g/3AT ss 10 cm/sec^ 
= 0.1 m/sec . The atmospheric Froude number is then on the order of 1.3. The Froude 
number exerts a strong influence on microburst shape while the Reynolds number affects 
turbulent details of the flow. In the laboratory, the Froude number can be modeled quite 
well, but the Reynolds number cannot. Let V^ A g = 100 cm/sec, g/JAT = 100 cm/sec^, 
v = 0.01 cm 4- /sec, and D = 0.5 cm. Then Fr = 14, and Re = 5000. Hence very high Froude 
numbers, or any smaller value can be achieved so that the shape of the microburst can 
be accurately simulated in the laboratory. Turbulence in the flow will not correspond 
quantitatively to that in the atmosphere but trends should be observable. 

Intuitively, the higher the Froude number, the shallower the outflow depth. 

The shallower outflow depth corresponds to an increased outflow velocity. The increased 
outflow velocity and decreased depth creates a zone of intense shear which results in 
increased turbulence production. Hence the Froude number has a secondary influence on 
turbulence. The Reynolds number affects turbulence and therefore turbulent mixing and 
thus the momentum diffusivity. Higher values of momentum diffusivity result in 
smoother velocity profiles (except near surfaces where viscous effects predominate) and 
thus the Reynolds number has a secondary effect on microburst shape. 

Figure 38 depicts a laboratory apparatus for studying microbursts. Saline solu- 
tions with densities 19% greater than fresh water can be achieved at 20°C. A reasonable 
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working value for 0AT would be 10%. For normal laboratory work, one meter head 
seems reasonable. This gives a maximum velocity of V=j2& = 4 m/sec. If the 
injector is constructed with a contraction, higher values of velocity can be achieved. The 
> valve controls the strength of the momentum source. The apparatus of Figure 38 can be 

of any size and complexity depending on resources. The valve could be microprocessor 
controlled with flowmeter feedback, for example. Velocity measurements could be made 
with hot wire or split film, laser Doppler velocimeter, or by particle photography 
methods. Flow visualization could be achieved by means by Schlieren, shadowgraph, 
interferometer, particle photography, hydrogen bubbles, or dye. 




t) it) 

Many of the questions posed previously could be studied with the apparatus of 
Figure 4. The effect of topography on microburst spreading should be modeled quite 
well with the apparatus since Froude number scaling can be achieved. Some idea of the 
t effect of surface roughness can be accomplished also. Surface roughness affects the 

- surface layer immediately but eventually could affect high., layers. Consider for the 

moment Figure 39. The laboratory microburst center is at the center of the 'ircle. 

The microburst front after a time is retarded by the increased roughness in the rough 
sector. The retardation creates shear layers in the transition area between the rough and 
smooth surfaces. If the effect is marked enough, Kelvin-Helmholz instabilities could 
be created. Variations of vorticity are observed in blowing dust associated with micro- 
hm ts. 



Figure 39. Effects of surface roughness on microburst spreading. 


By increasing head and hence velocities, variation of Reynold’s number can be 
achieved in the laboratory model. The effect on turbulence characteristics could be 
observed for changing Reynolds number. The turbulence variations, i.e., length scales 



and intensities have crucial relevance to the wind simulation model. While differences 
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will be observed between the laboratory turbulence and atmospheric turbulence, the 
trends should be there. The importance of the trends or functional forms was 
explained in Chapter V. 
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APPENDIX B. THE THREE-DIMENSIONAL CHARACTER OF THE 
VON KARMAN SPECTRA 


As an example, consider the three-dimensional von Karman spectrum for the 
vertical velocity component. 


d> 33 _ 


Cj (i^ 2 + i> 2 2 ) 


[1 + C 2 (y x 2 + v 2 2 + i' 3 2 ) 17 / 6 


(78) 


where 


440?r 


C, = -4— a*" A L 


2 A , 5 


and 


C 2 = (2naLY 


Defining p 2 = v j 2 + i^ 2 , substituting into (78) and solving for gives the following 
result. 


"3 = 1 
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i 2 
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(79) 


Equation (79) is the equation for surfaces of constant 033 . 

If Equation (79) is differentiated with respect to p, set to zero, and solved for p 
the result is 


/_ 6 _\ 17/22 / fl _\ 3 /’ 7 

\l 7 C 2 / U33/ 


( 80 ) 
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For this value of p for a given value of 4 > 33 , achieves its maximum value. 

The maximum value of < 1*33 can also be determined for given Cj and C 2 . By 
inspection of (78), the maximum value of 4*33 occurs at ^3 = 0, since f 3 is a moni- 
tonically decreasing function of $ 33 . The problem of finuing a maximum for <{>33 is 
reduced to finding the maximum value of $33 defined by 


4> 


C| 


33 ~[1 + c,p¥W> 


(81) 


Differentiating (81) with respect to p equating the result to zero, and solving for p gives 



Substituting this into the expression for ^33 gives 
a _ 6 Cj/l 1C 2 

Using 0 = 1 m/sec, L = 500 m in Equations (82) and (83) gives 

p = 1.76 x 10'^ nW and #33 = 1.37 x 10^ m^./sec^ 

The maximum value of 4*33 occurs at a small but nonzero value of p and has a large 


(83) 


but not infinite value. 
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APPENDIX C. SUM OF COMPLEX EXPONENTIAL SERIES 


In this study the following sum occurs frequently in derivations. 


N-l 

S = £ (e i0 ) n (84) 

n-"0 


The above summation is recognized immediately as a partial sum of the geometric series. 
Equation (84) can l e rewritten in the following form. 

oo oo oo oo 

s = £ (e i0 ) n - £ (e i0 ) n = £ (e i0 ) n - e i0N £ (e i0 ) n 
n=0 n=N n=0 n=0 

OO 

S=2> 0 )" [1 -e i0N ] (85) 

n=0 

Equation (85) is true if the infinite series converges. From complex variable theory it is 
known that the geometric series below convenes if |z| < 1. 


OO 



( 86 ) 


In actual fact (86) converges everywhere in the closed icgion |z| < 1 except at z * J: i . 
Hence 


S = 


1 - e 


i0N 


1 - e 


i0 


0 # 0, ff 


( 87 ) 
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By direct substitution into (84) S can be calculated at 0 = 0, and rr. 
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S = N 0=0 

S = 0 0 = 0, N even 
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APPENDIX D. DERIVATION OF PROPERTIES OF TRANSFORMED NOISE 


The purpose of this appendix is the derivation of the expressions summarized in 
Table 3. In these derivations, x n and x n 3 are assumed to be real white noise 


such that E[x n ] - 0, E[ x n ^] °x n “’ ^l' ( n x m^ ®mn> E ^ 


EJ x n 1 “,n 2 ,n 3 1 a x n “- and ^nj.n-i^m^m^iT^ a x n " 5 nimj 5 m 2 n 2 5 m 3 n 3 - 
The operator EM means the expected value of the variable in the parentheses and 5 mn 

is the Kronecker delta, i.e., 5 rnn = 0 for n, and 5 mn = 1 if m = n. 


Expression 1. 

E[X k ] =0 
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X k = S x n e~i 2 ™ k / M Ar 
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E[X k ] = E 
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M-l 

- £ E[x n ] e - j 27rnk / M = 0 
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Expression 2. 

M-l M-l 

E[X k X k -l = £ v E|x n| x„ 2 l e 
n l=0 n 2 =0 


■j27r(n r n 2 )k/M 2 
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Expression 3. 

E[X k X,*] = E 


Expression 4. 
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Expression 5. 


Efim 2 


[X k n “V* 
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m 2 Y 
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pression 6. 
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E(Re[X k lIm[Xjl 1 = Y Y E ( x ni x n? ] C0S < 2,rn l k / M > sin (2»« 2 1/M) (Ar) 2 
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nj=0 n 2=0 
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= a 2 (Ar) 2 y cos (2irnk/M) sin (2irnl/M) = 0 
n=0 

By the orthogonality property of the sine and cosine functions. 
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Expression 8 : 
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I exp 


1 (t!|-ti||k| (nj-njlki * n 3~ n 3> k 3 1{ 

“ m7~" m 7 ~ + m 3 Jj 


= f s ' (if| Ari Ar jf 


M , - 1 Mj-I Mj -1 


"i 


=0 ni-0 113=0 


6.1,7,, »n 2 *i 6 n 3 a 3 e**» 


| in|-n | )kj lnii!;)k" (nj-Kjjkj 

~vTj + Mi m 3 


; °x n ‘ (Arj Ari Ar 3 )- M, Mi M3 


Expression 9: 


E,X k|*2*3 X f, = °*n‘ (Af l Ar ' 


Mi-l Mi-' Mj-l j- n,(k.-l|) ti2 (k 2‘ 1 2 ) n 3 (k V , 3*"l 

Ar 3 >' III ~mT + mTJ 

npO n 2 =0 n 3 =0 L 


1 1 t s , 5. i M, M ; M 3 

= (Ar, Ar : top °x n ‘ 6 k,M ^2 3 


= [Mj M 2 M 3 *»52 V s 3 )2 ' * k 2 , 2 6k 3 1 3 


¥ 






, *4 


Vi* 
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Expression 10: 


Mi -1 M -,-1 


_ , 2itk->m Jirkjnj 

ElR e : tx lt , ik; . k ,)i=<Ar 1 £ 2-„ .5, tus ' <2tn ‘ k ' ,M ' — 


n^=0 n >=0 nj“0 


Mj-I M : -l M 3 -l 

= (Ar, Ar-, Ar 3 o r J] I I 
n l=0 n 2 =0 03=0 


1 \ [ / n l k l n 2 k 2 n ^ k A 

rr°'[ 4, \^ t V'sr/ 


= (Af| Ar> *ir 3 o x )" 


M, M : M, , 


T R C £ c *r 

n ] =0 


: (^rj <ir 2 Ar 3 o x r 


Mj M-, M 3 


1 i ^6 0k) + 6M r : (« 0 k : + 5 ~r k :^ ^ 8 ok, + 5M 3 - k .^ 


Expression 1 1 : 


M r l M 2 -l M3-I 

EI' m 2 (X k k , k )l = <Ar, AriAr^) 2 £ £ £ sin 2 

nj=0 m=0 n 3 = 0 


Mj -1 M 2 -l M 3 -l 

(Ar, Ar ; Ar 3 a^) 2 £ £ £ 

nj=0 n 2 =0 n 3 =0 


( n l k l n-,k 2 n 3 k 3 \ 

~m7 +- mT + "mJ J 

( n l k l n 2 k 2 ” 3 ^ 

1^ + -M7 + -Mij 


(Arj Ar 2 Ar 3 a x „)~ Mj Mi M 3 


(%,*\ (*•**«. j 


k 


* I ”0 • 






Expression 12 : 
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E l Re < x k,k 2 k 3 > 1 m< X l,Ul 3 )l = <°x n Ar l 4r 2 Ar 3> 2 


Mj-I M-,-1 

I t 


n l = ^ n 2 “^ 


(o Xn Ar, Ar 2 Ar$) : 
4 J 


Mj-| M-»-l 

Z Z 


nj^O n 2 a O 




+ 


exp 



n i(* 2 _k 2 ^ 

^ — — + 
Mi 


n 3 (l 3 -kj) 


- exp 



n 2 ^ 2 +, 2 i 

Mi 


♦ 


n 3 (k 3 +’ 3 )\ 

m, y 


(o Ar , ^r-. Ar 2 >- 

4 ] * k l>| **2': 5 k 3 l 3 + 4 0k, 5 0l, 4 0k 2 4 0b % 3 6 oi 3 

" 4 k, | i 4 k 2 l 2 4 k 3 ! 3 - 4 0k, S 0!| 4 0k 3 4 01 2 4 0k 3 4 01 3 
= 0 J 


k 


r-. & 
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APPENDIX E. ESTIMATE OF THE ONE-DIMENSIONAL SPECTRUM FROM 
THREE-DIMENSIONAL SIMULATED TURBULENCE 


This Appendix derives the result presented in Equation (65). The notation = 

X- Mj - 1 

Ai>j Ai >2 Av^ is used and the convention/ ,is accepted to meanX^ 


ki 


kj=0 


w(nj,n 2 ,n 3 ) = A 3 v EEE H(kj,k 2 ,k 3 ) N(kj,k2,k 3 ) 


h 1 
H M 

e L 


l k l n 2 k 2 
+ + 


M- 



k l k 2 k 3 


(89) 


In the above expression N(kj,k 2 ,k 3 ) is the transformed noise and H(kj,k 2 ,k 3 ) is the 
filter function corresponding to the desired output spectrum. 

Define W as follows. 


Qn , V' / ^ -j2ir(l 1 n 1 /N 1 ) 

W(lj,n 2 ,n 3 ) = 2^ w(ni,n 2 ,n 3 ) e 11 1 Arj 


(90) 


n l 


Substituting Equation (89) into Equation (90) and rearranging gives 


W = A 3 * Arj EEE H(k,,k 2 ,k 3 ) N(k 1 ,k2,k 3 ) e 

k, k 2 k 3 


f n 2 k 2 n 3 k 3l fniCkl- 1 !)" 

‘ [“2 + “sjy'e L N > ■ 


j2?r 


W = A 3 . Arj £ £ H(l 1 ,k 2 ,k 3 ) N(l 1 ,k 2 ,k 3 )e 
k 2 k 3 


f n 2 k 2 t n 3 k 3 1 

L M 2 + M 3 J 


With this result the variance of W becomes 


E[|W(l 1) n 2 >n 3 )| 2 ] =(A 3 j> Arj a^) 2 EE |H(lj,k 2 ,k 3 )| : 

k 2 k 3 



\ 

J 

t 

\ 

! 

f 

i 

i 


i 


(91) 


where a ^ is the variance of the transformed noise and is given in Table 3 as 



M 1 m 2 m 3 °n 2 


( ’»1 \ v 


( 92 ) 


The filter function is given in Equation (60) and is repeated here for convenience in 
digital form. 


H(lj »k 2 »k 3 ) 



v $2 % <1> w< 1 l> k 2’ k 3> 


a 


n 


(93) 


Here H(lj,k 2 ,k 3 ) and 4> w (lj,k2,k3) are taken to mean H(lj Ai^^Ai^^A^) and 
<b w (l| Ai>j ^At^l^A*^) respectively. 

Using Equations (92) and (93) in Equation (91) along with relations between sampl- 
ing frequencies and grid spacing yields the following desired result. 


E[|W(l 1)n 2,n 3 )l 2 l = 0 w (lj) = *w(ll*2<k 3 ) *>2 A "3 

1 k 2 k 3 


( 94 ) 



APPENDIX F. FROZEN TURBULENCE GENERATION IN THE 
FREQUENCY DOMAIN 




FTURB generates a three-dimensional block of turbulent velocities for either the 
u(II=l), v(II=2), or w(II=3) components. Dimensionless sampling frequencies are 
selected for each coordinate direction, = FS1, v s ^ - FS2, and v = FS3 in the pro- 
gram. A typical value for these parameters is 50. The first executable statement in 
FTURB is a call to a library routine LGBUF. The HP computer system this program was 
run on has a limited I/O buffer and a call to LGBUF creates a buffer LBUF 128 - 16 
bit words long so that records of this length can be read by the system. For other 
systems without buffer problems the call to LGBUF can be removed along with the 
dimensional array LBUF(128). 

The file of transformed turbulence is created in a direct access disc file in records 
1 28 - 1 6 bit words. The array being created is complex, and each complex number is 
made up of four words, so each block contains 32 complex numbers. Storage can be 
visualized as a one-dimensional array with index IX changing ;stest, 1Y, next fastest, and 
IZ slowest. Figure 40 depicts storage for a small block. In the block, each cell corres- 
sponds to a storage location (logically) for one complex number. The 128 words are 
stored physically in 4-32 complex word records 

Generation of the turbulence in the frequency domain requires access to lines of 
the block of logical turbulence. For example, in the figure, say we want to get a line of 
turbulence in the Z direction, say IX = 3, IY = 2, and IZ = 1 to 4. These indices corre- 
spond to words 1 1, 43, 75, and 107. Routine FCHRP handles this chore. 

Within the main program, the turbulence is generated according to the symmetry 
relations of Table 2. Variances of the transformed turbulence real and imaginary parts 
come from expressions 10 and 1 1 in Table 3. The (FI, 0,0) loop corresponds to 
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symmetry relation 5 in Table 3. The (0,F2,0) loop corresponds to relation 6, the 
(0,0, F3) loop to 7, the (F1.F2.0) loop to 8, the (F1.0.F3) to 9, the (0,F2,F3) loop to 
10, and the (F1,F2,F3) loop to expression 11. 


LOGICAL STORAGE 



RECORD 1 



PHYSICAL STORAGE 


Figure 40. Example of logical and physical storage for program FTURB. 
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S T«00004 IV OH CKOOOVV UV1HU 0««4* «U«« «>«««« 

FTN4X,L 
tFILESC 0,2) 

PROGRAM FTURB 

C 



C PROGRAM FTURB: THIS PROGRAM CREATES TRANSFORMED 3-D TURBULENCE BASED 
C ON THE VON KARMAN SPECTRAL MODEL. ITS PURPOSE IS TO CREATE THE FR02- 
C EN TURBULENCE DATA BASE FOR THE TURBULENCE SIMULATION MODEt 
C 

C VARIABLE DEFINITION: 

C Ml MAX - MAXIMUM NUMBER OF FI FREQUENCIES 

C N2MAX - MAXIMUM NUMBER OF F2 FREQUENCIES 

C N3MAX - MAXIMUM NUMBER OF F3 FREQUENCIES 

C FSf - FI DIRECTION SAMPLING RATE < DIMENSIONLESS > 

C FS2 ■ F2 DIRECTION SAMPLING RATE < DIMENSIONLESS > 

C FS3 - F3 DIRECTION SAMPLING RATE (DIMENSIONLESS) 

C ICART - LOGICAL UNIT NUMBER OF THE DISC DRIVE 

C II ■ VELOCITY COMPONENT C1-U, 2-V, 3-W > 

C 

C 

DIMENSION LABL< 10), LBUF( 1 28 ) 

COMPLEX X( 512 ), Y( 51 2 > 

COMMON N1H,N2M,N3M, INPUT, ICART 
INTEGERS N1 M , N2M , N3M 
C 

C THIS HP LIBRARY ROUTINE IS NECESSARY TO INCREASE I/O BUFFER SIZE. 

C 

CALL LGBUFCLBUF, 128) 

INPUT ■ 1 
ICART - 34 
ROOTS - SQRT< 2 . ) 

WRITE< INPUT, 9999) 

9999 FORMATC 1 1 HOUTPUT LU-? > 

READ< INPUT, 9998) LUOUT 
9998 FORHAT< 14 ) 

MRITE< INPUT, 9997) 

9997 FORMATC 7HN 1 NAX»? ) 

READC INPUT, 9998) N1MAX 
WRITEC INPUT, 9996) 

9996 FORMAT< 7HN2MAX-? ) 

READ< INPUT, 9998) N2MAX 
WRITE< INPUT, 9995) 

9995 F0RMAT< 7HN3MAX»? > 

READC INPUT, 9998) N3MAX 
WRITEC INPUT, 9994) 

9994 FORMATC 5HFS1 ■? ) 

READC INPUT, 9993) FS1 
9993 FORMATC FI 0.0) 

WRITEC INPUT, 9992) 

9992 FORMATC SHFS2-?) 

READC INPUT, 9993) F62 
WRITEC INPUT, 9991 ) 

9991 FORMATC 5HFS3-?) 

READC INPUT, 9993) FS3 
WRITEC INPUT, 9986) 

9936 FORMATC 38H ENTER VELOCITY COMPONENT Cl, 2, OR 3>> 

READC 1,9985) II 
9985 FORMATC I 1 > 
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FILE TO BE OPENED) 


NtH • HI MAX 
N2H • N2HAX 
N3M ■ N3MAX 

HREC « N1MAX*N2HAX*N3HAX/32 
HREC - N1M-N2H-N3IV32 
UR I TEC 1 ,8887) HREC 

6887 FORMAT < 7HHAXR2C- ,17, 25H ENTER 
REAOC t , 8806 ) < LABLC I ), 1 — 1, 10) 

6886 FORMATC 1 0A2) 

OPENC I CART , F I LE-LABL , I OSTAT- 1 08 , STATUS- ' NEW \ ERR-99 , RECL-256 , 
* FORM- 'UNFORMATTED ' , ACCESS- 'DIR ' , MAXREC-MREC ) 

DFf - FSl/’FLOATCNIMAX) 

DF2 - F82/FL0ATC N2MAX ) 

DF3 ■ FS3/FL0ATC N3MAX > 

N1D2 - N1MAX72 
N2D2 - N2MAX/2 
N3D2 - H3MAX/2 
DDDF - DFt*DF2*DF3 
XI MAX • FLOATC N1MAX) 

X2MAX - FLOATC N2MAX) 

X3MAX - FLOATC N3MAX) 

FFF - FSf-F$2-FS3 

SIGXK - SQRTCX1MAX*X2MAX*X3MAX/2. >/FFF 
P - PHIUC 0, , 0. , 0, , II ) 

SUM - P 


C 

C 

C 


H - SQRTCFFF-P) 

XC 1 ) - SIGXK*R0GT2*H*CMPLXC GRANC ),0, ) 
<F1 , 0/ 0 > LOOP 


DO 25 Kt - 2, N1 02 

FI - < FLOATC XI >-1 . )-DFI 
P » PHIUCF1 , 0 . , 0 , , 1 1 ) 

SUM - SUM +2 . -P 
H - 80RTC FFF-P > 

XCKO - SI GXK-H-CNPLXC GRANC ), GRANC >> 
XCN1MAX-K1+2) - C0NJGCXCX1 ') 

25 CONTINUE 

FI - FI ♦ DFI 
P - PHIUCF1 , 0. , G. , II ) 

SUM - SUM ♦ P 
H ■ 8QRTC FFF*P ) 

XCN1D2M) - SIGXK-R00T2-H-CMPLXC GRANC ), 0 . > 
CALL Ft tfPC2,t,NtMAX,l,l,1,t,X) 

C 

C C 0,F2, 0) LOOP 
C 


DO 30 K2 - 2,N202 

F2 - < FLOATC K2 )-1 . >-DF2 
P - PHIUC 0 . *F2, 0. , 1 1 > 

SUM - SUM ♦ 2 . *P 
H - 8QRTCFFF*P) 

XCK2) - SIGXK-H-CMPLXC GRANC ), GRANC >> 

XC N2MAX-K2+2) - CON JGC XC K2) ) 

30 CONTINUE 

F2 - F2 ♦ DF2 
P • PHIUC 0. ,F2,0. ,11) 

SUM - SUM ♦ P 
H - SQRTC FFF*P ) 

XCN2D2-M) - SIGXK*ROOT2*H*CMPLXCGRANC >,0. ) 


i 

j 

\ 


i 


t 


i 



i 

1 

i 

t 

i 


i 

i 






c 

C <0, 

c 


35 


C 

C <Ft 
C 


40 


45 


47 


C 

C < FI 
C 
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CALL FCHRPC2, 1 , t , 2,N2HAX, 1 , 1 ,X<2>> 

),F3) LOOP 

DO 35 K3-2,N3D2 

F3 » <FL0AT<K3>-1 . )<*DF3 
P - PHIU< 0. , 0. , F3, II > 

8UH - SUM ♦ 2.*P 
H - SQRT<FFF*P> 

X<K3> - SIGXK*H+CMPLX(GRAN< ), GRANC >> 

X< N3MAX-K3+2 > - C0NJCCX<K3>> 

CONTINUE 

F3 - F3 ♦ DF3 

P - PHIU< 0. , 0. ,F3, II > 

H • SQRT<FFF*P> 

X< N3D2+1 ) « SIGXK*ROOT2*H*CMPLXCGRANO, 0. > 

CALL FCHRPC2/ 1 , * t . I , 2, N7MAX, X< 2 ) > 

F2,0> LOOP 

DO 45 K1*2#N1D2 

FI - (FLOATCK! >-1 . >*DF1 
DO 40 K2 - 2, N2MAX 

F2 « <FL0AT<K2>-1 . >*DF2 
P - PHIUCF1 ,F2,0. , II ) 

SUM * SUM + 2.*P 
H * SQRT<FFF*P> 

X< K2 > - SIGXM-H^CMPLXCGRANOjGRANO) 

Y< N2MAX-K2+2 > - C0NJGCXCK2>> ' 

CONTINUE 

CALL FCHRP< 2 , K 1 , K 1 , 2 , N2MAX , 1 , I , X< 2 > ) 

CALL FCMRPC2,N1MAX-K1*2,N1HAX-Kt r2,2,N2MAX,1,1,Y<2>> > 

CONTINUE 
FI - FI ♦ DF1 
DO 47 K2-2,N2D2 

F2 - CFLOATC K2>-! . >*DF2 
P - PHIU<F1,F2/0. , II ) 

SUM - SUM ♦ 2 . *P ; i 

H = SQRT<FFF*P> 

X<K2> - SIGXK*H*CMPLX< GRANC >,GRAN< )) |T 

X< N2MAX-K2+2 > - C0NJGCXCK2>> 

CONTINUE > 

F2 ■ F2 ♦ DF2 1 

P * PHIUCF1 ,F2, 0. , II ) 

SUM - SUM ♦ P | 

H » SQRTCFFF*P> 

XCN2D2-HJ - SICXK*R00T2*H*CMPLX< GRANC 0 . ) 1 

CALL FCHRP<2,N1D2*1,N1D2*t,2,N2MAX,1,1,X<2>> l 

0,F3> LOOP | 

DO 55 K1*2,N1D2 j 

FI « CFLOATCK1 )-l . )*DF1 

DO 50 K3*2,N3MAX 1 

F3 » CFLOATC K3>-1 . >*DF3 I 

P - PHIUC FI , 0. , F3, II > 

SUM - SUM ♦ 2.*P 

H - SQRTCFFF*P> ! 

XCK3> « SIGXK*H*CMPLXC GRANC >, GRANC >> ! 

Y< N3MAX-K3+2 > ■ C0NJG<X<K3>> < 
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50 CONTINUE 

CALL FCHRP<2,K1 ,K1 , 1 , 1 ,2,N3HAX,X<2>> 

CALL FCHRP<2,N1NAX-KH-2,N1NAX-K1*2, 1,1 , 2, N3HAX, Y< 2 > > 

55 CONTINUE 

FI « FI ♦ DF1 
DO 57 K3 - 2,N3D2 

F3 - <FLOAT<K3>-1 . >*DF3 
P - PHIU<F1 , 0 . , F3, II ) 

SUN - SUN ♦ 2.*P 
H - SQRT<FFF*P> 

X< K3 > - SIGXK*H*CNPLX<CRAN< ),CRANO) 

X<N3NAX-K3+2> ■ C0NJG<X<K3>> 

57 CONTINUE 

F3 » F3 ♦ DF3 
P * PHIUtFI , 0. ,F3, II ) 

SUN - SUN ♦ P 
H » SGRT<FFF*P> 

X< N3D2+1 ) * SIGXK*R00T2*H'*CNPLX<GRAN< ),D.) 

CALL FCHRP<2,N1D2+1 , N1D2+1 , t , 1 , 2, N3HAX, X< 2 > ) 

C 

C <0,F2,F3> LOOP 
C 

DO 65 K2-2,N2D2 

F2 - <FL0AT<K2>-1 . )*DF2 
DO 60 K3»2, N3HAX 

F3 ■ <FL0AT<K3>-1 ■ >*DF3 
P » PHIU< 0. , F2,F3, II ) 

SUN - SUN ♦ 2.*P 
H • SQRT<FFF*P> 

X< K3 > « SIGXK*H*CNPLX<GRANO,GRAN< )) 

Y< N3NAX-K3+2 > - C0NJG<X<K3>> 

60 CONTINUE 

CALL FCHRP<2,1,1,K2,K2,2,N3HAX,X<2>> 

CALL FCHRP<2, 1, 1 ,N2HAX-K2+2, N2HAX-K2+2, 2, N3NAX, Y< 2 > > 

65 CONTINUE 

F2 - F2 ♦ DF2 
DO 67 K3-2,N3D2 

F3 = <FL0AT<K3>-1 . >*DF3 1'. 

P - PHIU<0.,F2,F3,II> - jf 

SUN - SUN ♦ 2.*P ” 

H - SQRT< FFF*P > 

X<K3) > SIGXK*H*CNPLX<GRANO,GRAN< >) , 

X< N3HAX-K3+2 > - C0NJG<X<K3>> 

67 CONTINUE j 

F3 =* F3 ♦ DF3 ! 

P - PHIU<0.,F2,F3,II> i 

SUN - SUN ♦ P * 

X< N3D2+1 > » SIGXK*H*R00T2*CNPLX<GRAN< ),0. ) ! 

CALL FCHRP< 2 , 1 , t , N2D2* 1 , N2D2+ 1,2, N3MAX , X< 2 > > ! 

C 

C < F 1 ,F2,F3 ) LOOP j 

C i 

DO 300 K1»2,N1D2 j 

FI » <FLOAT<I<1 >-1 . >*DF1 
DO 200 K2«2,N2NAX 

F2 - <FL0AT<K2)-1 . )*DF2 
DO 100 K3-2,N3NAX 

F3 ■ < FLOAT< K3 >-1 . >*DF3 
P - PHIU<F1 , F2,F3, I I > 

SUN - SUN + 2,*P 
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H ■ SQRT< FFF*P > 

X< K3 > - S1GXK*H*CHPLX<GRAN< ),GRAN< >> 

Y< N3MAX-K3+2 > ■ C0NJC<X<K3>> 

100 CONTINUE 

CALL FCHRP<2,X1 ,KI , K2, K2, 2, N3MAX, X< 2 ) > 

CALL FCHRP< 2, N1 MAX-K1 +2, Nt NAX-K1 *2, N2NAX-K2+2, N2MAX-K2+2, 

* 2,N3MAX,Y<2>> 

200 CONTINUE 
300 CONTINUE 

FI • FI ♦ OF1 
DO 350 K2-2.N2D2 

F2 ■ <FL0AT<K2>-1 . >*DF2 
DO 325 K3-2.N3MAX 

F3 - <FL0AT<I<3 >-1 . >*DF3 
P - PHIU<FI,F2.F3,II> 

SUN • SUN ♦ 2.*P 

X<K3> - SICXK*H*CMPLX<GRAN< >.GRAN< >> 

Y< N3MAX-K3+2 > ■ C0NJG<X<K3>> 

325 CONTINUE 

CALL FCHRP< 2, N1 D2+1 , N1D2+1 .K2,K2.2,N3NAX.X<2)> 

CALL FCHRP< 2.N1D2-H , N1D2+1 , N2MAX-K2+2 . N2MAX-K2+2 , 2 , N3NAX , Y< 2 > > 
350 CONTINUE 

F2 - F2 ♦ DF2 
DO 375 K3-2.N3D2 

F3 « <FL0AT<K3>-1 . >*DF3 
P - PHIU<F? .F2.F3. II > 

SUM - SUN + 2 , *P 
H - SQRT<FFF*P> 

X< K3 > • SICXK+H*CMPLX<CRAN< >.GRAN< >> 

X<N3NAX-K3+2> « C0NJG<X<K3>> 

375 CONTINUE 

F3 - F3 *• DF3 
P - PHIU< Ft , F2.F3, II) 

SUM * SUM ♦ P 
H • SQRT<FFF*P> 

X<N3D2+t > ■ SIGXK*H*CMPLX<GRAN< >.0. > 

CALL FCHRP< 2, N1 D2-M , N1 D2+1 , N2D2+1 . N2D2+1 .2. N3MAX, X< 2 > > 

C 

VAR » SUMoDDDF 

UR I TE< LUOUT , 9990 > N 1 MAX , N2MAX , N3MAX 
9990 FORMAT< 7H N1MAX-, I4.2X.6HN2MAX-, I4.2X.6HN3HAX-, I4> 

URITE< LUOUT. 9989 > 

9989 F0RMAT< / > 

WRITE< LUOUT. 9988 > FS1.FS2.FS3 

9988 FORMAT< 5H FS J = , F 1 0 . 4, 2X, 4HFS2*, FI 0 . 4, 2X, 4HFS3-, FI 0 . 4 > 

URITE< LUOUT, 9989 > 

WRITE< LUOUT, 9987 > VAR 
9987 FORMAK5H VAR-.F10.4) 

CLOSE< ICART, I0STAT-I0S.ERR-99, STATUS- 'KEEP ' > 

STOP 

99 WRITE< 1,8888 > I0S 

8888 FORMAT< "IOSTAT ERRORfi "I4> 

STOP 

END 

C 

C 

C FUNCTION PHIU: 

C THIS FUNCTION CALCULATES THE 3-D VON KARMAN SPECTRAL FUNCTION PHI 1 1 
C IN NONDIHENSIONAL FORM 
C 
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C F 1,F2,F3 ARE THE 3 COHPONENTS OF SPATIAL FREQUENCY 
C TLS IS THE TURBULEHT LEHGTH SCALE 

C 

C 

FUHCTION PHIU< FI , F2,F3, II > 

A « 1.339 
A4 ■ A**4 
PI a 4.*ATAH< 1 . > 

FSQ - FI*F1 F2*F2 + F3*F3 
ARC ■ FSG*< 2 . *PI*A )**2 
DENOM » < 1 . ♦ ARC >**< 1 7 . /6 . > 

IF< I I .EC. 1 > THEH 
FF » FI 
ELSE 

IF< I I .EC. 2) THEH 
FF ■ F2 
ELSE 

FF - F3 
END IF 
ENDIF 

PHIU - (440. /9. >*PI**3*A4*<FSQ-FF**2 VDENOM 
RETURH 
END 
C 

C 

C SUBROUTINE: FCHRP < FETCH REPLACE) THIS ROUTINE READS OR URITES 
C THE FOLLOWING < X< N1 , N2, N3 >, N1 =N1 MIN, N1 MAX ) 

C OR < X< N1 , N2, M3 >, N2=N2NIN, N2MAX ) 

>' C OR < X< HI , N2, N3 >, N3“N3MIN, N3MAX ) 

C FROM A RANDOM ACCESS DISC FILE. 

C 

1 c 

i C 

i SUBROUTINE FCHRP< IFOR, HI L, N1 H, N2L, N2H, N3L ,N3H, T > 

| COMMON N1N,N2M,N3M, INPUT, ICART 

COMPLEX T< 1 >,ZBUFF<32) 

INTEGERS IL, IH. N1 M, N2M, N3M 
' INTEGER DIDN 

C 

IL ■ N1L+<N2L-1 )*N1M+<N3L-1 )*N2M*N1M 
IH « N1H-KN2H-1 >*N1 M-K N3H-1 >*N2M*N1M 
C 

NRECL ■ < I L— 1 V32 ♦ 1 
NRECH - < IH-1 )/32 ♦ 1 
C 

IBL ■ IL - < NRECL- 1 >*32 
C 

DIDN * 0 



IFCN1L .NE. N1H > 

DIDN - 

1 


IF< N2L .NE. N2H > 

DIDN = 

N1M 


IF< N3L .NE. N3N) 

DIDN - 

N1M*N2M 


IF< DIDN .HE. 0) GO TO 25 



WR1TE( INPUT, 9999) 


9999 

FORMAT* 13HERROR 

FETCH 

1 > 

c 

STOP 



25 

IS ■ 1 




18 - IBL 

NSTEP • DIDN,' '32 

IF< NSTEP .LT. 1> NSTEP « 1 


I 



i 

5 

i 

♦ 

i 

I 

i 


l 
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00 too IREC - NRECL , NRECH , NSTEP 

READ< ICART,REC-IREC, I0STAT-I0S,ERR-99> < ZBUFF< 1 >, I-t ,32> 
50 IF< IFOR .EC. t> T<IS> « ZBUFF< IB > 

IF< IFOR .EC. 2 ) ZBUFF< IB > - T<IS> 

IS ■ IS ♦ t 

IB * IB ♦ DION 

IF< IB ,LE. J2> GO TO 50 

IB « N0D<IB,32) 

IF< IB .EO. 0> IB - 32 

IF< IFOR .EQ. 2) WRITE< I CART, REC* IREC, I0STAT**I0S,ERR“99> 

* <ZBUFF< I ), 1*1 ,32> 

IflO CONTINUE 
RETURN 

99 MR I TE< INPUT, 9998 > I OS 

9998 F0RHAT< 19HFETCH IOSTAT ERROR ,I3> 

STOP 

END 

END# 
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APPENDIX G. TRANSFORMATION OF FREQUENCY DOMAIN 
“TURBULENCE” TO THE SPACE DOMAIN 


Program FFT3D transforms the file created by FTURB to the space domain via 
the inverse FFT. FFT3D contains a call to library routine LGBUF which was explained 
in Appendix F. The transform is calculated by operating on small lines (one-dimensional) 
from the data in the previously created file. This program is extremely memory efficient 
but time inefficient. Core storage is minimal but numerous disc file reads and writes 
are required. 

The three-dimensional FFT is calculated using a simple one-dimensional FFT 
routine. The three-dimensional inverse FFT is defined by 


(n i k i n 2 k 2 n 3* c 3\ 

M, M 2 M 3 / 

- 1 


(95) 


This equation is broken up into three groups of one-dimensional transforms as follows 


m 3 -i 


1 j27T n^ki/M? 

A(k l5 k 2 >k 3 )=— 2- X(ki ) k 2 ,k 3 ) e 33 3 

1 k 3 =0 


(96) 


The above procedure must be performed Mj x M 2 times and is implemented in the K3 
loop of the main program. The next operation is given in the following equation 


B(k 1 ,n 2 ,n 3 ) = ^ 


m 2 -i 

1 £A» lW . i! " 21 * 

n 2 =° 


(97) 


This transform is performed Mj x M 3 times. The third operation is given by 
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x ( n l>n 2 ,n 3 ) = — 


1 


M r l 

L 


kj =0 


B(kj,n 2 ,n3) 


}2ir njkj/Mj 
c 


(98) 


and is performed M 2 x M3 times in K1 loop. For a 64 x 64 x 64 array only 64 complex 
words must be in memory at a given time. 
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3D T-00004 18 OH CR00033 USING 00837 BLKS R-0000 

FTN4X,L 
8FILES< 0,2) 

PROGRAM FFT3D 

C 

C 

C PROGRAM FFT3D i THIS PROGRAM CALCULATES A THREE DIMENSIONAL FFT IN 
C PLACE FOR DATA RESIDING IN MASS STORAGE < TYPE 1 DIRECT ACCESS FILE ) 

C 

C 

DIMENSION LABL< t 0),L8UF< 126) 

COMPLEX X< 31 2 ) 

COMMON N1HAX,N2MAX,N3MAX, INPUT, ICART,K1 ,K2,K3 
INTEGERS N1 MAX, N2NAX, N3MAX 
CALL LGBUF<LBUF, 128) 

INPUT - 1 
ICART - 34 
R00T2 ■ SQRT( 2 . ) 

WRITE< INPUT, 9999 > 

9999 FORMAK tIHOUTPUT LU-?) 

REAO< INPUT, 9998 > LUOUT 
9998 FORMAT< 14 > 

URITE< INPUT, 9997) 

9997 FORMAK 7HN 1 MAX*? > 

RE AD< INPUT, 9996) N1MAX 
URITE< INPUT, 9996) 

9996 FORMAK 7HN2MAX-? > 

READ< INPUT, 9998) N2MAX 
WRITE< INPUT, 9995) 

9995 FORMAK 7HN3MAX-? ) 

READ< INPUT, 9998) N3MAX 
WRITE< INPUT, 9994 > 

9994 FORMAK SHFS 1-?) 

READ< INPUT, 9993) FS1 
9993 FORMAK FI 0.0) 

URITE( INPUT, 9992) 

9992 FORMAK 5HFS2-?) 

READ< INPUT, 9993) FS2 
URITEC INPUT, 9991 ) 

9991 FORMAK 5HFS3-?) 

READ< INPUT, 9993) FS3 
MREC - N1 MAX*N2MAX*N3NAX/32 
WRITE< 1 ,8887) MREC 

8887 FORMAK 7HMAXREC-, 17, 25H ENTER FILE TO BE OPENED) 

READ< 1,8886) <LABL< I ), 1-1 , 1 0) 

8886 FORMAK 10A2) 

OPEN< I CART , F ILE-LABL , I OSTAT- I OS , STATUS- ' OLD ' , ERR-99 , RECL-256 , 

* FORM* 'UNFORMATTED ' , ACCESS- 'DIR ' , MAXREC-MREC > 

XI MAX - FLOAKN1MAX) 

X2MAX • FLOAK N2MAX ) 

X3MAX - FLOAT< N3MAX) 

DF1 - FS1/X1MAX 
DF2 - FS2/X2MAX 
DF3 - FS3/X3MAX 
N1D2 - N1 MAX/2 
N2D2 - N2MAX/2 
N3D2 - N3MAX/2 
DDDF - DF1*DF2+DF3 
FFF - FS1 *FS2*FS3 
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1 

/■ 


HIM - NfHAX 
N2M « N2MAX 
N3M » N3MAX 
WRITE* 1,6999) 

6999 FORMAT* -BEGIN K3 LOOP") 

C 

C K3 LOOP 
C 

DO 200 K 1 ■ 1 , N 1 MAX 
DO 100 K2*1 , N2MAX 
WRITE< 1,6998) K1 ,K2 

6998 FORMAT* " K1-"I3" K2«"I3 H |A"> 

CALL FCHRP* 1 , XI , K 1 , K2, K2, 1 , N3N, X > 
CALL FFT*X,N3MAX, 1 > 

CALL FCHRP*2,K1,K1,K2,K2, t,N3M,X> 
100 CONTINUE 
200 CONTINUE 
C 

C K2 LOOP 
C 

URITE< 1.699?) 

6997 FORMAT* "BEGIN K2 LOOP") 

DO 400 K1-1.N1MAX 
DO 300 K3-1.N3HAX 
WRITE* 1,6996) K1 ,K3 

6996 F0RMAT*"K1«“I3" K3»"I3"tA"> 

CALL FCHRP< 1 , K 1 . K1 . 1 , N2M, K3, K3, X > 
CALL FFT<X,N2MAX, 1 ) 

CALL FCHRP<2,K1 ,K1 , 1 , N2M, K3, K3, X ) 
300 CONTINUE 
400 CONTINUE 
C 

C XI LOOP 
C 

WRITE* 1,6995) 

6995 FORMAT* "BEGIN K1 LOOP") 

DO 600 K2»1 , N2HAX 
DO 500 K3-1.N3HAX 
URITE< 1,6994) K2.K3 

6994 FORMAT* "K2-" 13“ K3-"I3"^A"> 

CALL FCHRP* 1 , 1 ,N1M,K2,K2,K3,K3,X) 
CALL FFT< X,N1MAX, 1 ) 

CALL FCHRP< 2,1 ,N1M,K2,K2,K3,K3,X) 
500 CONTINUE 
600 CONTINUE 
C 

C MULTIPLY BY CONSTANT 

C 

WRITE* 1,6993) 

6993 FORMAT* “MULTIPLY BY CONSTANT") 

CONST - FFF 
WRITE* 1,5554) CONST 
5554 FORMAT* “C0NST«6"E1 2 .5 ) 

DO 950 K3-1.N3NAX 
DO 900 K2-1.N2MAX 

CALL FCHRP* 1 , 1 , N1 M, K2, K2, K3, K3, X ) 
DO 800 J«1 , N1MAX 
X* J ) ■ X* J)*CONST 
800 CONTINUE 

CALL FCHRP* 2, 1 ,N1M,K2,K2,K3,K3,X) 


i 



i 


i 

I 

i 


! 


i 


I 
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900 CONTINUE 
950 CONTINUE 

CL0S£< I CART , STATUS- * KEEP ' , I OSTAT- I OS , ERR-99 ) 

STOP 

99 WRITE< 1,5555 > IOS 

5555 FORMAT< 17HI0STAT ERROR NO. ,14) 

STOP 

END 

C 

C — 

C SUBROUTINE: PCHRP (FETCH REPLACE) THIS ROUTINE READS OR WRITES 
C THE FOLLOWING < X< N t , N2, N3 >, N1 -1 , N1 NAX ) 

C OR < X< N 1 , N2 , N3 > , N2- 1 , N2HAN > 

C OR < X< N1 , N2, N3 >, N3-1 , N3MAX > 

C FROM A RANDOM ACCESS DISC FILE. 

C 

C 

C 

SUBROUTINE FCHRP< IFOR , N l L , N1 H, N2L , H2H , N3L , N3H, T > 

COMMON N 1 MAX , N2MAX , N3MAX , INPUT, ICART,K1 ,K2,K3 
INTEGER*4 N 1 MAX , N2MAX , N3MAX , IL, IH 
COMPLEX T< I ), ZBUFF< 32 > 

INTEGER DIDN 
C 

IL - N1L-KN2L-I )-N1MAX*<N3L-1 )*N2MAX*N1MAX 
IH - N1H-KN2H-1 )*N1NAX«-<N3H-1 >*N2HAX*N1 MAX 
C 


C 

C 


9999 

C 

25 


50 


500 


NRECL - (IL-I)/^ ♦ 1 
NRECH - (IH-1V32 ♦ 1 


IBL - IL - < NRECL-1 )*32 


DIDN - 0 

IF(N1 L .HE. N1H) DIDN « 1 
IF< N2L .HE. N2H) DIDN - N1MAX 
IF< N3L .NE. N3H > DIDN - N1NAX-N2MAX 
IF< DIDN .NE. 0) GO TO 25 
URITE< INPUT, 9999) 

FORMAT< 13HERR0R FETCH 1) 

STOP 


IS - t 

IB - IBL 

NSTEP ■ DIDN/32 

IF< NSTEP .LT. 1 ) NSTEP » 1 

DO 100 IREC • NRECL, NRECH, NSTEP 

READ< ICART,REC-IREC, I0STAT-10S, ERR-99) < 2BUFF< I >, 1-1,32) 
IF< IFOR .EQ. 1) T< IS ) - ZBUFF(IB) 

IF< IFOR .EG. 2) ZBUFF(IB) - T<IS) 

IF< IB . GT . 32 ) GO TO 500 

IF< IB.LT. 1 ) GO TO 500 

IF< IS.LT. 1 > GO TO 500 

IF( IS . GT .512) GO TO 500 

IF<CABS<T< IS >>.GT. 1 .E5> GO TO 500 

IF<CAB3<ZBUFF< IB)). GT . 1 .E5> GO TO 500 

GO TO 600 

WRITE< 1,7998) K1,K2,K3 
MRITE< 6, 7998 > K1,K2,K3 
WRITE< 1 , 7999 > IB, ZBUFF< IB >, IS, T< IS) 

WRITE< 6, 7999 ) IB,ZBUFF< IB), IS,T< IS) 




1 

1 

• 

! 

* 

I 

I 

l 

i 


l 
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7999 FORNATC ZBUFF< "15“ >**2<X,£I2.5>* T< "15" >-"2< X,EI2,5>> 

7998 FORMATC " K1--I3* K2-“13 H K3-"I5> 

600 IS » IS ♦ I 

IB - 18 + DIDN 

1F< IB .LE. 32 > CO TO 30 

IB - MOCX IB, 32) 

IF< IB .EC. 0) IB - 32 

IF< IFOR .EC, 2) WRITE< ICART , REC- IREC , IOSTAT- IOS , ERR-99 > 

* <ZBUFF<I),I«1,32> 

tOO CONTINUE 
RETURN 

99 WRITE< INPUT, 9998) IOS 

9998 FORMAT< 19HFETCH IOSTAT ERROR ,I3> 

STOP 

END 

C 

C 

C SUBROUTINE) FFT 

C JIM COOLEY'S SIMPLE FFT PROGRAM— USES DECIMATION IN TIME ALGORITHM 
C X IS AN N-2*<*H POINT COMPLEX ARRAY f HAT INITIALLY CONTAINS THE INPUT 
C AND ON OUTPUT CONTAINS THE TRANSFORM 

C THE PARAMETER INV SPECIFIED DIRECT TRANSFORM IF 0 AND INVERSE IF 1 

C 

C 

SUBROUTINE FFT< X, N, INV) 

COMPLEX X< 1 > 

COMPLEX U, W, T, CMPLX 
INTEGERS N 
C 

C X - COMPLEX ARRAY OF SIZE N— ON INPUT X CONTAINS 

C THE SECUENCE TO 8E TRANSFORMED 

C ON OUTPUT X CONTAINS THE DFT OF THE INPUT 

C N - SIZE OF THE FFT TO BE COMPUTED— N«2*+M FOR 1.LE.M.LE.I5 

C INV - PARAMETER TO DETERMINE WHETHER TO DO A DIRECT TRANSFORM <INV«0> 

C OR AN INVERSE TRANSFORM < INV-O 

C 

M - ALOC< FLOAT< N ) )/ALOG< 2 . ) ♦ .1 
NV2 - N/2 
NM1 ■ N - 1 
J - 1 

DO 40 I-I,NM1 

IF U.GE.J) CO TO 10 
T « X<J> 

X<J> - X<I> 

X<I) ■ T 
10 K • NV2 

20 IF <K.CE.J> GO TO 30 

J - J - K 
K - K/2 
GO TO 20 

30 J ■ J ♦ K 

40 CONTINUE 

PI . 4 . *ATAN< 1 . > 

DO 70 L»t ,M 
LE - 2**L 
LEI » LEZ2 
U ■ <1 .0,0.0) 

U ■ CMPLX<COS<PIZFLOAT<LE1 > >, -SIN< PIZFLOAT< LEI >>) 

IF < INV.NE. 0) W ■ CON JG< W ) 

DO 60 J»1 ,LE1 
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do 30 i«j,n,le OF POOR QUALITY 

ip - i ♦ LEI 
T - X<IP>*U 

x<ip> - x<n - t 



X< I > - X< I > ♦ T 

i 

50 

CONTINUE 

* 


U - U*U 


60 

CONTINUE 


70 

CONTINUE 

1F< INV.Eft.O) RETURN 
DO 60 I-t,N 

X<I> • X< I >/CMPLX< FLORT< N >, 0 , > 


60 

CONTINUE 

RETURN 



END 

END* 

) 
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APPENDIX H. FLIGHT OF A POINT “AIRPLANE” THROUGH 
WINDSHEAR AND TURBULENCE 




& 


Program JAWS2 moves a point through the JAWS wind shear and generated 
turbulence at a constant selectable ground speed. The “airplane” can fly a straight 
horizontal, three degree glide slope approach, or three degree glide slope departure f ■ 
path. Initial positions within both the JAWS data and the turbulence are selectable 
through flight mode variables MODEF. Two real variable “knobs”, JCON and TCON 
permit an increase in turbulence or turning on or off of the wind shear so that separate 
components can be examined separately. 

Interpolation procedures for this program were described in Chapter VII. The 
interpolation for the JAWS data is handled by subroutine NTERP. Functions for the 
calculation of turbulent length scales and gust intensity are handled by functions TLS 
and SIGX, respectively. The demonstration functional forms of TLS and SIGX are also 
described in Chapter VII. 

Storage for the JAWS data is in cell form as described in Chapter VII. The 
turbulence is stored in 16 bit integers and are converted to real numbers by dividing 
by 10000. Storage of the JAWS data is on disc unit 34 and the turbulence on disc 
unit 37. 



j 
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2 T-00004 ie ON CR00033 USING 000SS 8LK8 R-0000 



FTN4X,L 
SFILES< 0>2> 

PROGRAM JAW82 


PROGRAM JAUS2i THIS PROGRAM CALCULATES WINDS FOR A PLANE 
MOVING THROUGH THE JAWS DATA *ET. 


PROGRAMMER) WARREN CAMPBELL 


VARIABLE DEFINITION: 

NBLPF • NUMBER OF BLOCKS PER FIELD 

NBLPL - HUMBER OF BLOCKS PER LINE 

2 LINE - LINE NUMBER OF CURRENT CELL 
IBL1K - PRESENT BLOCK NUMBER 

NX1 ■ NUMBER OF CELLS IN THE X DIRECTION 

NYI ■ NUMBER OF CELLS IN THE V DIRECTION 

NZt - NUMBER OF CELLS IN THE Z DIRECTION 

XMIN - MINIMUM VALUE OF X <KM> 

XMAX ■ MAXIMUM VALUE OF X <KM> 

YMIN » MINIMUM VALUE OF Y <KM> 

YMAX » MAXIMUM VALUE OF Y <KM> 

ZMIN - MINIMUM VALUE OF Z <M> 

ZMAX - MAXIMUM VALUE OF Z <M> 

DELX - X GRID SPACING <M) 

DELY - Y GRID SPACING <M> 

DELZ - Z GRID SPACING <M> 

DELT - TIME STEP CSEO 

IX - X INDEX OF CURRENT CELL 

IY - Y INDEX OF CURRENT CELL 

IZ » Z INDEX OF CURRENT CELL 

I XT - X TURBULENCE INDEX 

IYT • Y TURBULENCE INDEX 

IZT - 2 TURBULENCE INDEX 

ITPT ■ NUMBER OF CURRENT TURBULENCE POINT 

XT » DIMENSIONLESS TURBULENCE X LOCATION 

YT ■ DIMENSIONLESS TURBULENCE Y LOCATION 

ZT » DIMENSIONLESS TURBULENCE Z LOCATION 

N1MAX - MAX NUMBER OF TURBULENCE POINTS IN X DIRECTION 

N2MAX - MAX NUMBER OF TURBULENCE POINTS IN Y DIRECTION 

N3MAX > MAX NUMBER OF TURBULENCE POINTS IN Z DIRECTION 

FS1 - DIMENSIONLESS SAMPLING FREQUENCY IN FI DIRECTION 

FS2 » DIMENSIONLESS SAMPLING FREQUENCY IN F2 DIRECTION 

FS3 - DIMENSIONLESS SAMPLING FREQUENCY IN F3 DIRECTION 

TCON » TURBULENCE CONSTANT. LARGE VALUE ELIMINATES TURBULENCE. 

JCON - JAWS CONSTANT. 1. FOR WIND SHEAR ♦ TURBULENCE 

0. FOR TURBULENCE ONLY 

X - CURRENT X LOCATION <KM> 

Y « CURRENT Y LOCATION <KM) 

Z - CURREN1 2 LOCATION <M> 

VX - EAST-WEST INERTIAL VELOCITY <P0SIT1VE EAST IN M SEC> 

VY • NORTH-SOUTH INERTIAL VELOCITY < POSITIVE NORTH IN M SEC> 
VZ - VERTICAL VELOCITY (POSITIVE UP IN M/SEC) 

XO - INITIAL X VALUE < KM ) 

YO - INITIAL Y VALUE <KM> 

20 • INITIAL Z VALUE (M> 

HDG ■ HEADING < DEGREES > 

VCRS - GROUND SPEED (H/SEC) 
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VAIRS - AIRSPEED (M/SEC) OF POOR QUALITY 

UWINO - EAST-WEST WIND SPEED (POSITIVE EAST IN H/SEC) 

/WIND » NORTH-SOUTH WIND SPEED (POSITIVE NORTH IN M/SEC) 
WWIND - VERTICAL WIND SPEED (POSITIVE UP IN NHSIC> 

MODEF * FLIGHT NODE 

1 • HORIZONTAL FLIGHT 

2 • ILS APPROACH (3 DEGREE GLIDE SLOPE > 

3 - TAKEOFF (3 DEGREE GLIOE SLOPE > 

CELL ■ ARRAY CONTAINING CORNER POINTS FOR PRESENT FIELD IN 
PRESENT CELL 

I CELL - CURRENT CELL NUMBER 

ICWBL = CURRENT CELL NUMBER WITHIN BLOCK (1-16) 

ICB * BEG1NING INDEX OF CELL WITHIN BLOCK (1,9, 17, ETC. ) 

ICE * ENDING VALUE OF CELL WITHIN BLOCK - ICB+7 
XC - X VALUE WITHIN CELL (0-DELX> 

YC - Y VALUE WITHIN CELL ( 0-DELY > 

ZC - Z VALUE WITHIN CELL ( 0-DELZ ) 

T » TIME 

NPTS * MAXIMUM NUMBER OF POINTS TO BE CALCULATED. THIS IS A 
MEANS OF TERMINATING THE PROGRAM EARLY FOR DEBUGGING 
PURPOSES. 

IPT ■ NUMBER OF THE PRESENT POINT BEING CALCULATED 
SIGMX * MAXIMUM VALUE OF GUST INTENSITY 


DIMENSION LBUFF( 128), IA( 1 28 >, NAHE( 1 0 >, LA8ELC 40 >, CELL( 8 >, ITT( 128> 
COMMON DELX,DELY,DELZ,D3, SIGMX 
REAL JCON 

INTEGERS ICELL, ITPT,N1MAX,N2MAX,N3HAX 
CALL LGBUF( LBUFF , 1 28 > 

PI « 4 . *ATAN( 1 . ) 

WRITE( 1 , 9999 > 

9999 FORMAK 'ENTER FILE TO BE OPENED*") 

RE AD( 1,9998) < NAME< I >, I“1 ,10) 

9998 FORMAK 4 0A2) 

WRITE( 1 , 9997 > 

9997 F0RMAT( "ENTER FLIGHT MODE*"/ 

* " 1 - HORIZONTAL FLIGHT"/ 

* * 2 - ILS APPROACH (3 DEGREE GLIDE SLOPE)"/ 

- " 3 * TAKEOFF (3 DEGREE GLIDE SLOPE)") 

READC 1,*) MODEF 
WR I TEC 1,9996) 

9996 FORMAK "ENTER CONSTANT GROUND SPEED*") 

READC 1,*) VGRS 
WRITEC 1,9995) 

9995 FORMAK "ENTER TRUE HEADING IN DEGREES*") 

READC 1 ,*) HDC 
WRITEC 1 ,9994) 

9994 FORMAK "ENTER INITIAL POSITION (X0,Y0,Z0) IN (KM,KM,M)*" ) 

READ( 1 ,*> X0, Y0,Z0 
WRITEC 1,9988) 

9988 FORMAK "ENTER DELTA T IN SECONDS*") 

READC 1,*> DELT 
WRITEC 1,9993) 

9993 FORMAK "ENTER TAPE HEADER (40A2)*") 

READC 1 ,9998) ( LABELC I >, I«1 , 40 > 

WRITEC 1,9979) 

9979 FORMAK "ENTER MAXIMUM NUMBER OF POINTS*") 

READC 1,*> NPTS 
WRITEC 1,8999) 

8999 FORMAK "ENTER JCON,TCON*"> 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 
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8998 

8997 

8995 
7777 

9992 

9991 

9990 

9989 

8994 

9937 

8996 


9986 

9985 

9984 

9983 


READ<1,*> JCON, TCON 
WRITE* 1,8998 > 

T ORMAT* "ENTER N 1 MAX , N2MAX , N3MAX£ " > 

READ* 1 , • > N1 MAX, N2MAX, N3MAX 
WRITE* 1,8997) 

FORMAT* "ENTER FS1 , FS2, FS3£* ) 

READ* 1 , * ) FS i , FS2, FS3 
XTMAX » N1NAX/TS1 
YTMAX = N2MAX7FS2 
ZTMAX * N3MAX/FS3 
WRITE* 1,8995) 

FORMAT* "ENTER XTO, YTO, ZT0$" ) 

READ* 1 , * > XTO, YTO, 2T0 
WRITE* 1,7777 > 

FORMAT* "ENTER SIGMXA" ) 

READ* 1 , * > SICMX 

WRITE* 8, 9998 > < LABEL* I ), 1*1 ,40) 

WRITE* 8, 9992 > MODEF 
FORMAT* • FLIGHT MODE - "I2> 

WRITE* 8, 9991 ) HDG 

FORMAT* “ TRUE HEADING * “F6.1" DEGREES") 

WRITE* 8, 9990) VGRS 
FORMAT** INERTIAL HORIZONTAL VELOCITY-“F8. 0" M/SEC" > 

WRITE* 8,9989 > X0,Y0,Z0 

FORMAT*" INITIAL POSITION! X0="F8.2" kK Y0="F8.2" KM Z0=“ 
* F6. 0* M") 


WRITE* 8, 8994 ) XT0,YT0,2T0 

FORMAT*" XT0-"F8 .3“ YT0-"F8.3" ZT0="F8.3> 

WRITE* 8, 9987) DELT, TCON 

FORMAT** TIME INCREMENT - “F8.3" SECONDS TC0N = "E12.4> 
OPEN* 34, FILE-NAME, STATUS- 'OLD ' , IOSTAT-IOS, ERR-99 
■* , ACCESS- 'DIR ' , RECL-256, MAXREC-7201 ) 

WRITE* 1,8996) 

FORMAT* "ENTER TURBULENCE FILE NAMEfi" ) 

READ* 1 ,9998) < NAME* I ), 1-1 , I 0 > 

OPEN* 37, FILE-NAME, STATUS- 'OLD ' , IOSTAT-IOS, ERR-99 
* , ACCESS- 'DIR', RECL-256, MAXREC-2048) 

RFAD* 34, REC-1 , IOSTAT-IOS, ERR-99 > < IA< I >, 1-1 , 16 > 

NX1 * IA< 1 ) - 1 
NY1 - I A* 2 ) - 1 
NZ1 - IA<3) - 1 
XMIN - IA< 7 VI 00 , 

XMAX - IA<8)71 00. 

DELX - IA< 9 ) 

YMIN * IA*10V100. 

YMAX ■ IA< 11 >7100. 

DELY - I A* 12) 

ZMIN « I A* 13) 

ZMAX - I A* 14) 

DELZ • I A* 15 > 

D3 - DELX*DELY*DELZ 
WRITE* 8,9986 > XMIN, XMAX 
FORMAT*" XMIN-"F8 .2" KM 
WRITE* 8,9985 > YMIN, YMAX 
FORMAT* " YMIN-"F8 . 2" KM 
WRITE* 8,9984) ZMIN, ZMAX 
FORMAT* " ZHIN-"F8 . 2" KM 
WRITE* 8,9983) 

FORMAT*" DELX 
WRITE* 8,9960 ) 


XMAX*"F8 . 2" KM") 


YMAX-"F8 . 2" KM") 


ZMAX--F8.2" KM") 


DELX, DELY, DELZ 

"F6.0" M DELY-"F6 . 0" 

NPTS 


DELZ-"F6 , 0" M" ) 
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9960 FORMAU" NPTS-“I5> 

URITE< 8,9982 > 

9982 F0RHAT<8X“T , '9X"X“9X" Y"9X"Z"?X“UttIND“5X“VWIND“5X*yUIND*5X 

*“VAIRS N > 

q*** ******** CALCULATE • OF BLKS PER LINE AND • BLKS PER FIELD 
IF< HOD( NX 1 ,16) . EQ . 0) THEN 
NBLPL - NX 121 6 
ELSE 

NBLPL - NX 1 2 1 6 ♦ 1 
END IF 

NBLPF » NBLPL*NY1*NZJ 

C******** CALCULATE VX,VY,VZ ********* 

THETA * PI*HDG2180 . 

VX » VGRS*SIN< THETA) 

VY - VCRS*COS< THETA ) 

IF<H0DEF . EQ. 1) THEN 
VZ » 0. 

ELSE 

IF<N0DEF .EQ. 2) THEN 

VZ - -VGRS*TAN< 3 . *PI2180 . ) 

ELSE 

VZ * VGRS*TAN< 3 . *PI21 80 . ) 

END IF 
END IF 

C****** INITIALIZE POSITION AND TINE ******** 

X - XO 
Y - YO 
Z - ZO 
T - 0. 

IPT * 1 
XT * XTO 
YT » YTO 
ZT * ZTO 


C****** CALCULATE CURRENT CELL INDICES, IX, IY, IZ, IXT, ETC. ***** 
50 IX « 1 000. *<X - XMIN >2DELX ♦ I 


IY - 1000. *<Y - YMIN )2DELY + 1 
IZ » <Z - ZMIN>2DELZ ♦ 1 
IF< XT .GT. XTMAX > XT = XT - XTMAX 
IF< YT .GT. YTHAX) YT - YT - YTMAX 
IF< ZT .GT. ZTNAX > ZT - ZT - ZTHAX 

IF<XT .LT. 0.) XT * XT ♦ XTMAX 

IF< YT .LT. 0.) YT ■ YT ♦ YTMAX 

IF< ZT .LT. 0.) ZT - ZT ♦ ZTMAX 

IXT * IFIX< XT*FS1 ) ♦ 1 
IYT » IFIX< YT*FS2 ) + 1 
IZT ■ IFIX< ZT*FS3 > + 1 
D WRITE< 6, 5000 ) IXT, IYT, IZT 

5000 FORMAK" IXT-"I5" IYT«“I5" IZT«*I5> 

D WRITER 6, 5001 ) XT,YT,ZT 

5001 FORMATC " XT--F8.3* YT-"F6.3" ZT-"F8.3> 

ITPT - IXT ♦ < IVT-1 >*N1HAX ♦ < IZT-1 >*N2MAX*N3MAX 
!F<HOD< ITPT, 128) .EQ. 0) THEN 
ITBLK ■ ITPT2128 
ELSE 

ITBLK - ITPT21 28 ♦ 1 
END IF 
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I PUB • ITPT - < ITBLK-1 )*128 
D URIT£<6, 5605) ITPT, ITBLK 
5005 FORHATC • ITPT-“I8* ITBLK-M5) 

C***** CALCULATE UUIND **** 

ICELL - IX ♦ NX1*<IY-t> + <IZ-1 >*NXt*NY1 
ILINE « IY ♦ < IZ-I >*NY1 
IF< HOD< IX, t < > .EG. 0> THEN 
ICUB * 16 

IBLOK - < ILINE-1 )*NBLPL ♦ IX/M6 ♦ 1 
ELSE 

I CUB ■ MOD< IX, 16) 

IBLOK « <ILINE-1 )*NBLPL + IX/16 * 2 
ENDIF 

READ< 34, REC-IBLOK, IOSTAT«IOS,ERR«99) < IA< I >, 1-1 , 128) 
ICB - <ICU8-1>*8 ♦ 1 
ICE - ICB + 7 
IC * 1 

DO 100 IB-ICB,ICE 

CELL< IC > - IA< IB >/1 00 . 

IC - IC ♦ 1 
100 CONTINUE 

READ<37,REC-ITBLK, I0STAT-I0S,ERR-99> < ITT< I >, 1-1 , 128> 
IF< JCON .LT. 0.001 > THEN 

UTRB - ITT<IPUBV< 10000. *TCON > 

ELSE 

UTRB - ITT< IPUB)*SIGX<X,Y,Z V< 10000. *TCON) 

ENDIF 

XC « 1 000 . *< X-XMIN > - < IX-1 >*OELX 
YC « 1 000. *< Y-YMIN ) - <IY-1)*OELY 
ZC * Z - ZHIN - <IZ-1)*DELZ 
CALL NTERP< CELL, XC, YC, ZC, UUIHD > 

UUIND - JCON*UUIND ♦ UTRB 
D URITE<6,5003> SIGX<X, Y,Z),UTRB 
5003 FORHATC" SIGX»"F8.3" UTRB«"F8.3> 

C***** CALCULATE VUIND ***** 

IBLOK - IBLOK + NBLPF 

READC 34, REC-IBLOK, I0STAT-I0S,ERR-99) < IA< I ), 1-1,128) 
IC - 1 

DO 125 IB-ICB, ICE 

CELLCIC) « IA< IB )/1 00 . 

IC - IC ♦ 1 
125 CONTINUE 

CALL NTERP< CELL ,XC,YC,ZC, VUIND) 

O**** CALCULATE U ***** 

IBLOK - IBLOK ♦ NBLPF 

REAO< 34, REC-IBLOK, IOSTAT-10S, ERR-99) < IA< I ), 1-1 , 1 28) 
IC - 1 

DO 150 IB-ICB, ICE 

CELL< IC ) - IA( IB VI 00 . 

IC - IC 1 
150 CONTINUE 

CALL NTERP< CELL, XC,YC,ZC, UUIND) 

C***** CALCULATE AIRSPEED ***** 

VAIRS - VGRS - UWIND*SIN< THETA) - VUIND*COS< THETA > 


C****« URITE RESULTS ***** 
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WRITEt 8/9981 ) T, X, Y, Z , UWIHD, VWIHD, WUtND, VAIRS 
9981 FORMATt 5X,F6 . 2, 7t 2X, F8 . 2 > > 

C***** INCREMENT TIME AND POSITION ***** 

T - T «• DELT 
X • X ♦ VX*DELT71 000 . 

Y - Y ♦ VY*DELT/’1 000 • 

Z - 2 ♦ VZ*DELT 

XT - XT ♦ VX*DELT/’TLS< X, Y, Z > 

YT - YT ♦ VY*DELT/*TLSt X, Y, 2 > 

ZT - ZT + VZ*DELT/TLStX,Y,Z> 

IF< IPT .LT. NPTS > THEN 
IPT - IPT ♦ 1 
GO TO 250 
ELSE 

ENDFILE 8 

CLOSEt 34, STATUS- 'KEEP ' , I OSTAT-IOS, ERR-99 > 

CLOSE< 37, STATUS- 'KEEP ' , IOSTAT-ICS, ERR-99 > 

STOP 

CNDIF 

250 IF< < X . GT . XMIN .AND. X.LT.XMAX) .AND. 

* < Y . GT . YMIN .AND. Y.LT.YMAX) .AND. 

* (Z.GT.ZMIN .AND. Z.LT.ZMAX>> THEN 

GO TO 50 

ELSE 

ENDFILE 8 

CLOSEt 34, STATUS- 'KEEP ' , IOSTAT-IOS, ERR-99 > 

CLOSEt 37 , STATUS- ' KEEP ', I OSTAT- 1 OS . ERR-99 > 

ENDIF 

STOP 

99 URITEt 1,9980) IOS 

9980 FORNATt "IOSTAT ERROR ifi"I4> 

END 

FUNCTION TLSt X, Y, Z ') 

C 

C FUNCTION TLSt TURBULENT LENGTH SCALE FUNCTION. FUNCTIONAL VALUE 
C IS FROM ED41 TERRESTRIAL ENVIRONMENT DOCUMENT. 

C 

TLS - 31 .5*< 2/^18. 3>**0. 64 

RETURN 

END 

FUNCTION SIGXt X, Y, Z > 

C 

C FUNCTION SIGXt THIS ROUTINE CALCULATES THE U TURBULENCE 
C STANDARD DEVIATION FOR THE JAWS JULY 14,1982 CASE. 

C 

COMMON DELX,DELY,DELZ,D3,SIGMX 


C ********************************************* ************ 

C t XCTR, YCTR > ARE THE COORDINATES OF THE MICROBURST CENTER 
C IN KILOMETERS. 

C********************************************************* 
XCTR * 14.2 


YCTR - -1 .5 

R - SQRTt < X-XCTR >**2 - <Y-YCTR>**2> 
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C**************************************************************** 

C CALCULATE THE XY FACTOR. THIS IS BASED ON B-57B NEASUREHNTS 
C FROM THE SAHE DATE. THE 5? ENCOUNTERED THE UPPER PART OF A 
C HICROBURST AND EXPERIENCED A DECREASE IN TURBULENCE AS IT PASSED 
C THROUGH THE DOWNDRAFT, FOR THE SANE RUN THE AVERAGE OF SIGX,SIGY, 

C AND SIGZ WAS ABOUT 6 M/SEC . FXY CAUSES A DECREASE IN SIGX OF 
C ABOUT 6 AWAY FROH THE OUTFLOW CENTER, AND 3 IN THE CENTER. ANOTHER 
C TERN IN FXY CAUSES SIGX TO DECREASE TO EXP<-.5> AT A DISTANCE OF 
C 5 KN. THIS IS IN KEEPING WITH A DECREASE IN TURBULENCE AWAY FRON 
C WIND SHEAR. 

C*****«**4^*******************************************+************ 

FXY » < SIGNX - .5*SIGNX*EXP<-R*R7< .75**2*2. > > >*EXP< -R*R/50 . > 
IF < Z . GT . 300. > THEN 
FZ • EXP<< 2-300. )/200 . > 

ELSE 

FZ « 1. 

ENDIF 

SIGX - FXY*FZ 

RETURN 

END 

SUBROUTINE NTERP< CELL , XC , YC , ZC , V > 


INTERPOLATION ROUTINE 


DIHENSION CELLO > 

COhHON OELX, DEL Y,DEL2,D3, SIGNX 

PH1 1 - < DELX-XC )*< DEL Y-YC >*< DEL2-ZC >/D3 

PH 1 2 ■ XC*< DELY-YC >■*< DELZ-ZC )/D3 

PHI3 « XC*YC*< DEL2-ZC VD3 

PH 1 4 - (DELX-XC >*YC<*< DELZ-ZC >/D3 

PHIS - ( DELX-XC )*< DELY-YC >*ZC/’D3 

PH 1 6 - XC*< DELY-YC >#ZC/D3 

PH I 7 ■ XC#YC*ZC7D3 

PH 1 8 - < DELX-XC >*YC*ZC7D3 

V • PH1 1 *CELL< 1 ) ♦ PHI2*CELL< 2 > ♦ PHI3*CELL<3> ♦ PHI4*CELL<4> 
• + PHI5+CELL<3> ♦ PHI6*CELL< 6 > * PHI7+CELL<7> ♦ PHI8*CELL<8> 
RETURN 
END 
END* 



